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Introduction 


onse Of an intinitely long uniform circular cylindrical shell 
‘is investigated. The shell is immersed in an acoustie fluid 
ied to move in the dilatational, translational, and inexten- 
il modes only. The shock wave is assumed to have its wave 
» the axis of the shell. 
ilvsis for the circular evlindrical shell has been formulated 
using a modal method of analysis in which the deflexion is 
he sum of modes and each mode is then evaluated indepen 
method has been applied only to the first two modes of a 
fh is subjected to a plane shock wave of constant pressure, 
tuse the evaluation of the detlexion for each mode is tedious 
excessive when several modes are required. Mindlin and 
ve simplitied the analysis by an approximation in which each 
element of the shell is assumed to act as an element of a flat plate. This 
plane wave approximation is valid only for small values of the time, 


generally much less than the transit time which the wave takes to cross the 


evlinder. Although the approximation appreciably reduces the analysis 


+] 
the 


results are not accurate at the interesting times when the modal 
deflexions reach a maximum, occurring for some modes at times larger 
than the transit time. 

In this paper an approximate relation is developed between the fluid 
pressure and particle velocity of a evlindrical wave. Using this relation 
the modal equations of motion are simplitied and in the typical examples 
chosen lead to numerical results in close agreement with the true elastic 
response. The numerical calculations required for this cylindrical wave 
approximation are no more extensive than those required for the plane 
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SUMMARY 


The modal method of analysis is used to determine the response of an elastic 
cylindrical shell subjected to a pressure pulse propagating through the surrounding 
acoustic medium. An approximate relation is established between the fluid pressure 
and velocity of a cylindrical wave, which is then used to obtain an approximate 
solution for the elastic response of the shell. 


1. Introduction 

THE elastic response of an infinitely long uniform circular cylindrical shell 
to a shock wave is investigated. The shell is immersed in an acoustic fluid 
and is constrained to move in the dilatational, translational, and inexten- 
sional-flexural modes only. The shock wave is assumed to have its wave 
front parallel to the axis of the shell. 

The elastic analysis for the circular cylindrical shell has been formulated 
by Carrier (1) using a modal method of analysis in which the deflexion is 
expanded as the sum of modes and each mode is then evaluated indepen- 
dently. The method has been applied only to the first two modes of a 
cylinder which is subjected to a plane shock wave of constant pressure, 
primarily because the evaluation of the deflexion for each mode is tedious 
and becomes excessive when several modes are required. Mindlin and 
Bleich (2) have simplified the analysis by an approximation in which each 
element of the shell is assumed to act as an element of a flat plate. This 
plane wave approximation is valid only for small values of the time, 
generally much less than the transit time which the wave takes to cross the 
cylinder. Although the approximation appreciably reduces the analysis 
the results are not accurate at the interesting times when the modal 
deflexions reach a maximum, occurring for some modes at times larger 
than the transit time. 

In this paper an approximate relation is developed between the fluid 
pressure and particle velocity of a cylindrical wave. Using this relation 
the modal equations of motion are simplified and in the typical examples 
chosen lead to numerical results in close agreement with the true elastic 
response. The numerical calculations required for this cylindrical wave 
approximation are no more extensive than those required for the plane 


(Quart. Journ. Mech, and Applied Math., Vol. XI, Pt. 2, 1958] 
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wave approximation and for times larger than the transit time involve only 
tabulated Bessel functions. 


2. Equations governing the interaction between the shell and fluid 
A plane wave is considered travelling in the direction of the axis of x at 
right angles to the axis of the cylinder as shown in Fig. 1. 


Incident 
Radiated pressure ee 
pressure pulse Cylindrical 
pulse 





The radially inwards displacement w and transverse displacement v when 
confined to the dilatational, translational, and inextensional-flexural modes 
can be expressed in generalized coordinates in the formt 


w 
w= > q,,cosné, 


0 
> In sin nO. (2.1) 
ann 


The radiated fluid pressure P. and radial velocity U, (in the radially out- 
wards direction) created by the presence of the shell, are defined by the 


velocity potential ad as _ ed 


P. —_ rr’ , 
ae , ér 


(2.2) 
where the velocity potential satisfies the hydrodynamic equation 
A2 
cya = SF (2.3) 
et? 


+t Baron (3) has removed this assumption and has presented a method using the plane 
wave approximation which allows for extensional effects in all modes. These effects are 
found to have a damped oscillatory form and are of importance only very shortly after 
the arrival of the shock wave. 
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c being the velocity of sound, p the density of the fluid. The velocity poten- 
tial, pressure, and velocity are expanded in the form 


d => cos n6 


Pp. = ¥ Pp, cosnd U, = > U,,008n6 





P.=3P,,cosond U,=> Un 008 00 | 
n=0 n=0 

P, being the pressure of the incident shock wave and U; being the radially 

inwards component of the particle velocity in the incident shock wave at 

the surface of the shell. Using equations (2.1), (2.2), and (2.4) the equations 

of motion of an elastic shell take the form 


a2 
a ES 


op 
My “et? OR In a int Pp 2 


at 
= — Fn 
or’ 


atr=a (2.5) 


ee (1—n?)? EAS 

»~ 72(1—p2)a* 

and a is the radius, h is the thickness, m is the mass per unit area of the 
shell, Z is Young’s modulus, and y is Poisson’s ratio. 

The condition that the radial velocity is continuous at the shell is that 


Od, U, Opp, 
TS 


for n > 1, 


atr=a. (2.6) 
ot 


The generalized velocity potential ¢,, must satisfy the equation 
1 6 [ad n® 1 a 
= onli at hance Gk. ait am ements 2.7 
= ae) a Pn c* at* (2.7) 


The solution is obtained by solving (2.7) for ¢, and eliminating ¢,, from 
equations (2.5) and (2.6). By using a Laplace transformation the result 
can be expressed in contour integral form as 


P,,,(z)K;,(z)—peU ,,(z)K,,(2) 
af (A, 22-+0,,)K,(2)—2K,,(2) jer S 





(2.8) 


where K,, is the modified Bessel function of the second kind 


a = 


a pa’ 


pt » — mn dw? dK, 
: dz 


as Coca pc? ’ K;,(2) ” 
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and P,,,(z), U;,(z) are the Laplace transforms of P,,,, U;,, with respect to 7’. 
The numerical solution of equation (2.8) in each mode for given values of 
A,,, ¢, and incident shock wave involves an extensive computing programme. 


3. Relation between particle velocity and pressure 

3.1. Résumé of the plane wave approximation 

A first approximation to the solution of equations (2.5), (2.6), and (2.7) 
has been given by Mindlin and Bleich (2). In that approximation the solu- 
tion of equation (2.7) is replaced by the relation 


ap, ~— l Od, 3.1 
= —- 8, (3.1) 
or c ot 
which is equivalent to the approximation K‘,(z) = —X,,(z) in equation 
(2.8). The approximation thus assumes that each element of the shell sub- 
jected tothe shock wave radiates a plane wave into the surrounding medium; 
an assumption which is increasingly in error as the interaction proceeds. 


Using this approximation the equations of motion in (2.5) reduce to 


* n T € 
Gn + pe + wd, = Pin tpcUin- (3.2) 
* dt dt 
Mindlin and Bleich (3) give the solution of (3.2) for the case of a plane 
shock wave of constant pressure and they evaluate the first three modes 
for a typical shell. 


3.2. The cylindrical wave approximation 
In order to obtain a better approximation we return to equation (2.7). A 
solution of this equation is given by Lamb (4) in the formT 


4, = A, n(- =) do; do = [ flet—reosh w) du, 
0 


where A,, is an arbitrary constant. It is shown in the Appendix that an 
alternative solution of similar form is 


¢, = { F,(ct—r cosh u) cosh nu du, (3.3) 
0 


where the function F, is such that the integral is convergent. A straight- 
forward calculation using (3.3) leads to the relation 


@ 


; , 
On = _%__! | g,, F,,(ct—r cosh u) cosh nu du, (3.4) 
or cot 6 


0 


+ A diverging cylindrical wave is taken here. The subsequent analysis can be readily 
extended to a converging cylindrical wave or a combination of the two waves. 
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__ 1+4nsinhutanh nu 
a 1+coshu : 


Equation (3.4) immediately shows the relationship between the fluid 
pressure and velocity, reducing to the plane wave approximation of (3.1) 
when the integral term is zero, i.e. for very large values of r. 

In order to examine the physical significance of the integral in (3.4) and 
to obtain a better approximation than in (3.1) we consider the case n = 0. 
Introducing the transformation rcosh u = ,/(r?+-z*) in (3.3) and (3.4) we 
find 





where 


i 4) 


_ f Rfict—(r?+2%)} 
do = | 2 J@P+24) dz, 





(3.3 a) 


9) 


Obo _ _ > _ 1 f tn Fate as, 


or cot or 

22\)-1 
where %Q= ( +,/( +5)| 3 
If z is taken as the coordinate along the axis of the shell then Lamb has 
shown that ¢, is the resulting potentiai of a system of point sources distri- 
buted over the axis of z with uniform line density F(t), each elemental 
contribution being taken at the retarded time and decaying inversely with 
distance. In the higher modes line densities of a more complex nature can 
be similarly obtained. 

For cylindrical waves the integral term which appears in (3.4a) can be 
regarded as a line density g,(z/r)F,(t) taken at the retarded time and decaying 
inversely with distance. The line density changes along the axis as given 
by go(z/r) decreasing steadily from a value of $ at z = 0 to zero as z/r + 00. 

An approximation is now made that g,,(z/r) can be replaced by a constant 
j,- The integral in equation (3.4) then reduces to the integral for the 
velocity potential and the relation becomes 

On = On Ing or Un 1 oF, 


m= —2 In p , 3.5 
ér cat 6 et pe ot Vr oe atte 





(3.4a) 


Integrating from a time ¢, to a variable time ¢ at constant distance r we find 


t 
=eolB— salt) fo | P,,(t) dt. (3.6) 
to 


U,,(t)— Un (to) _ 


pe p 


If the time ¢, = 0 is chosen as the time of arrival of the disturbance it is 
seen that the velocity in the fluid at a later time is a function not only of 
the instantaneous pressure but of all previous changes in pressure after the 
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arrival of the disturbance. It follows that in a cylindrical wave the fluid 
is left with a final velocity, or after-flow, following the passage of the 
pressure wave, and this after-flow disappears only after the pressure falls 
below its equilibrium value. 

In the cylindrical symmetric diverging wave we have 


t 
U;(t)— Ut) = OF) , 90 [ Piey ae 
pc pr 
to 


and the after-flow term is just g, times the corresponding after-flow term 
in the well-known case of a spherical wave. 


4. Calculation of the after-flow coefficient 


It is noticed that the range of g,, is 


wa 
= oe 
<n forn>1 


where the limits correspond to the times t = 0 and too. The variation 
of g,, with z/r is shown in Fig. 2. 

As a rough value for g, the value of g,, at z = 0 could be used. It follows 
from Fig. 2 that for 7, = } a larger weight is always attached to the contri- 
bution from the line density at increasing distance z along the axis. For 
the higher modes a smaller weight is attached to the contribution from 
these elements. The after-flow term in (3.4) will therefore be increasingly 
in error for later times when the contribution from the line density of 
distant elements becomes important. In the present approximation there- 
fore a mean value for g, is taken. Using the transformation z = rtan @ and 
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integrating (g,,/7)tan@ from — 4 to $7 the mean values for the first few 
modes are 


§,= 0-727  gG, = 2-099 
Jo = 1-192 gj, = 2-545 


Thus, the mean value g, increases with increasing values of n. Using the 
above value for g, less weight is attached to the effect from elements of the 
line density at small distances along the axis and greater weight is attributed 
to elements at large distances along the axis. The reverse is true for the 
higher modes. 


5. The approximate equations of motion 


Having obtained an approximate relation between the U,,, and P.,, it is 
possible to eliminate these quantities from equations (2.3) and (2.4). The 
equation of motion of the shell for the nth mode becomes 


aP, dU sn) 
a\q7 +9 Pin + pe aT |’ 


13q - d*q dq, a 

r, @ Fn 4 (y 1) % 4a n 

n dT? ( n Int 1)- dT? +a n dT + on In Tn pc? 
(5.1) 
y ty 

where . = —, o, = —. 


pa pe® 


The solution of (5.1) which replaces the exact solution (2.8) is 





2 f [Fat 2ME+Gu)+zpeUin(2)\ por gy 


: Sipe An(2—24)(2—24)(2—2) 


where z,, z., and z, are the roots of the eqration 

A, B+(AaGratle+e,2+0,9, = 90. 
The numerical evaluation of equation (5.2) in contrast to (2.8) is straight- 
forward. It merely involves calculating the contribution from the poles 


at z,, Z,, and z,, and at any poles which may be contained in the numerator . 
of the integrand. 


6. A comparison of the approximations 

6.1. The suddenly applied hydrostatic pressure 

When a hydrostatic pressure is suddenly applied to the shell at time 
t = 0 only radial displacements occur and the incident particle velocity is 
zero. For thin shells the mass has a negligible effect on the response except 
for a very small time in the initial stages of the motion. Taking the mass 
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as zero, the solution of equation (5.1) for a suddenly applied pressure P, is 
w= I 
= fof) tetor| sin }7'(4g,.04—02)—c08 $7 (4G %— “| |. 
w5 | V(49o%—9%) j 
This response is shown in Fig. 3 for o, = 1-08 and g, = 0-363. The exact 
solution calculated from equation (2.8) and that obtained using the plane 


Exact solution—_l,_ — 
a 

/Approximate-, Z | 

solution y . 


; 


as ? —-—+-+ = + +- ——--4 
Solution using plane wave 
approximation 








{ 2 3 


ct 
a 


Fic. 3. Cylinder response for a suddenly applied hydrostatic pressure 


wave approximation are also shown. The comparison illustrates how the 
cylindrical wave approximation follows the exact response to within a few 
per cent at any instant, whereas the plane wave approximation cannot 
simulate the overshoot of deflexion above the static deflexion P,/w?. 


6.2. The incident exponential shock wave 
An exponential shock wave of the form 
P, = Pye-Pt2+90 ft (2+-a)/c} 

is taken, where f is the decay constant and H(t) is the Heaviside unit 
function. The transformed pressure and radial velocity with respect to T 
on n PB, e-] ) U __ En PB, es I’ (z) 
+B n ’ a pe 2+B. n\“/> 
where /,,(z) is the modified Bessel function of the first kind and 


pa (2) = Fal) _ {1 for n re 

c dz \2 forn > 0. 
The mass of the shell is of importance in the dilatational and flexural modes 
only for the early accelerated motion of the shell plating shortly after the 


arrival of the shock wave. Therefore the evaluation of the response at all 


P,,,(z) = 


in 


€ 
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times apart from the very early stages of the motion is appreciably simplified 
by neglecting the mass of the shell. The calculation of the response at very 
early times is obtained by an alternative procedure in a later section. 


91°2 Shock wave of 41+2 i 
constant pressure poe ok nag a 
—— shock wave of 
410 Exact Approximate 41: constant pressure 
solution solution 


NG ° 
Exponential SI 


06 shock wave 70° 


= © cec-! 
B= 754 8e°: P ; 
704 f Approximate solution 
for exponential shock 


wave 
4024 b =< -1 
4 B=75q 5° 











ae a oe 
et ct 
a a 


Fic. 4a. Dilatational mode Fic. 46. Flexural mode 








In the dilatational and flexural modes the solution of (5.1) neglecting the 
mass of the shell is 
_ -q m (Fn eee B)L,e(— B)— BI,,_16(— B)\ 
En ln Pf % \ (z,+B)(z,4+- B) J 
ie (Gn —n-+-2,)L,9(%)+% se e-"(—T) 4. 
* (z,+ B)(z,—2.) 
{(J,—+-22)1,9(22) +22 aed aaa for n ~ 1) 
LF BY, —2,) 
cos~*(1—T) 


1 cos nO e°9.d@ (for T' < 2), 
where [,o(z) = (7 
J (6.1) 


eBU-T) 4 








| 
a os 
! 





I,(z) (for T' > 2), 
and z,, z, are the roots of 
22+o,z+0, 9, = 0. 

Taking as typical parameters the values g, = 0-363, a) = 1-08, and 
J, = 1-19, o, = 0-0047, the dilatational and first flexural deflexions are 
shown in Figs. 4 (a) and 4 (6) respectively, for a shock wave of constant 
pressure and a shock wave having an exponential decay constant 

= © sgec-1, 
15a 
The exact solution calculated from (2.8) is also given for the dilatational 
mode. 
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In the translational mode n = 1 the stiffness of the shell plating w, is 
zero. The translational motion is controlled at all times by the total mass 
of the shell, in general the shell plating itself being only a fraction of the 


Shock wave 
i i 1-2 of constant 
Approximate solution rentmate 
for shock wave of. oe pressure 


constant pressure 1-0 
S—___ Exact 
solution 
08 


Exponential 
shock wave 
Approximate = 3 -1 
solution for B sae 
exponential 
shock wave 
=< -1 
B sa 28° 








4 








c 
a 


Fic. 5a. Translational mode Fic. 5b. Translational velocity 
total mass. A typical case is a shell which when immersed in an acoustic 
fluid is in a neutrally buoyant condition; we then have A, = 1. 


For n 1 the translational deflexion is 


pe” , ((91- - U)Jo9(0)- - 21,(0) “(+1 Yoo) + 
J 


ap," | 12,B 
(9,—1)(2 Byho() , of + B—di dhol — B)— Blog < =m T) 4 





+2=-——_— —— +2 


1aB ' (z+ B)B 


of (1 +91— Whol@1) +21 Lool41)\ --2a- (6,2) 
\ A; (2, + B)zi 


where z, (A, 9, +-1)/A,. For a neutrally buoyant shell A, = 1, the trans- 





i 
+ 


lational displacements are shown in Fig. 5 (a) for a shock wave of constant 
pressure and a shock wave having an exponential decay 8 = (c/15a) sec". 
The translational velocities are shown in Fig. 5 (b). The exact velocity 
distribution obtained by Murray (5) for a shock wave of constant pressure 
is also shown. 


7. Properties of the approximate response 
Equation (5.1) can be expressed in the form 


m, M,, d°q d*q dq 
a n ny l M mo 9 | h + ae 2 
pe dé (m,, a) dt? ' (m,, L,,)%p dt + Wh, Tn 


M,, dP in. yf diy 


= Fin dt dt 
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_ op M,, wn 

where M,, 5.’ a, = Spo(m, + M,)" 

The first term, depending on the rate of acceleration, is important only 
in the early stages of motion. The inertial term contains an equivalent 
mass ©, arising from the water flow in the vicinity of the shell. The third 
term is a damping term having the damping factor a,. M,, is inversely 
proportional to g, and therefore decreases with increasing modes. Since 
w® is proportional to (n?—1)? the damping factor increases rapidly with 
increasing » and modal deflexions become highly damped in the larger 
modes. Furthermore, the deflexion in the inextensional or flexural modes 
either oscillates about or tends exponentially to the static modal deflexion 
P,,,/w2 at large times. The condition for the occurrence of oscillations is 


2 2pc\* 
wt < (-) (m,+M,) (n #1). 
n 
For thin shells m,, can be neglected and the condition becomes 
o, < 49, (n #1). 
Taking og = 1-08 in the inextensional mode n = 0, oscillations occur if g, 
is greater than 0-27. This result is illustrated in Fig. 3 where the cylindrical 
wave approximation is calculated for g, = 0-363 and where the plane wave 
approximation corresponds to g, = 0; oscillations occur only in the cylin- 


drical wave approximation. For the flexural modes n > 2, the condition 
becomes 


(n*—1)" <6 (ita) jn (n> 2). 


Taking as a typical valuea/h = 100 and noting that g,, lies between } and n, 
oscillations disappear in some mode between 8 and 19 and in all higher 
modes. 


The limiting value of modal velocity and displacement for large times 
can be obtained from (2.8) without resorting to the cylindrical wave approxi- 
mation. In particular, the limiting translational velocity is 


dq “al : (8 = 0) 
im ( ‘Z) =v, = {pce\A,+1 
$i 0 (B #0). 
In the latter case the final displacement when the shell comes to rest is 


cals 
=|, 77 # 0). 
The approximate non-zero limiting values which are obtained from (5.2) are 
= By att) - 2 (att) 
AGi+1)’  peB\Ayg,+1 


lim 9g; = qo = 
t+oa 
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The approximate and exact values coincide for the neutrally buoyant 
cylinder A, = 1. They also agree when the value of 7, is taken equal to the 
asymptotic value of g, for large values of z/r, i.e. when g, = 1. In the latter 
case the correct weight is attributed to elements of the line source at large 
distances along the axis which because of retarded time effects can only 
influence the motion at large times. 

After the transit time 7' = 2 the translational velocity contains no oscil- 
latory terms and the motion tends to v,, as the sum of two exponential terms. 
For the neutrally buoyant cylinder the velocity which is derived from (6.2) is 


pe dq, a (22,(B) | 21,( B) \ pa-7)_ 24a(\21)) pina T>2 
Py dt | B " (\z,|— B)| (z,|—B) (1 > 2), 
where B = aB/c, z, = —(g,+1). 


For small values of B, B < |z,|, the second exponential obviously decays 
more rapidly than the first and therefore the translational velocity tends 
to zero from positive values. In the special case B =- 0, the velocity tends 
to v,, from below this value and the velocity at all times is less than v,,. 
The corresponding asymptotic behaviour of the exact solution is discussed 
by Murray (5). In particular when B = 0 the exact solution has a maximum 
velocity greater than v,, occurring at a finite time 7’~ 5. The above 
difference in the asymptotic behaviour is a consequence of the cylindrical 
approximation. 


8. The response at early times 

In order to determine the modal deflexions for values of the time appre- 
ciably smaller than the transit time, equation (5.1) must be solved without 
neglecting the mass of the shell. This solution will involve the evaluation 
of the incomplete Bessel function J,,9(z) for large negative values of z and 
for small times numerical inaccuracies become apparent. An alternative 
procedure is to develop a power series in time 7’. This is accomplished 
using the exact solution (2.8) by expanding the integral as a series in z~! 
and integrating each term. The leading terms of this series are 


og, = 8 (GT fap, 2 48] GTP | 
€,0, A te A,,Vr | P(3) ii 4 ri) + Me 


This expression is used to determine the deflexions for values of the time 
appreciably smaller than the transit time, 2a/c. 
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APPENDIX 
The solution of 4, given by (3.3) may be verified by substitution in the differential 
equation (2.7). We find 
c "bn . py nbn _ é *bn 


or?" r Or r? ect? 


- | { @F,(x) coshu OF, (x) _ Y 
= sinh? - _ -- 

\ ex 
ry 











; - . oz Ff, “5 Fa(z)}eosh nu du 
a r ex 








| cos 2 > 

[s ec miawen {é “Fn _ntF,| hes [<= oF, tM sinhnu Fy] 

| u® j a ou u=0 
0 


2 -osh nu sinh u é OF Az) 


r ox 
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where x = ct—rcoshu. The result follows sae F,(t) = 0 for eae values 
of ¢t exceeding a certain limit. 
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SUMMARY 

General formulae are obtained for the second-order effects in the deformation of 
incompressible transversely isotropic elastic bodies. 

The problem of the second-order torsion and extension of a homogeneous incom- 
pressible cylinder transversely isotropic with respect to its generators is reduced to 
the solution of the classical torsion and flexure boundary-value problems together 
with another boundary-value problem involving two complex potential functions. 
Without solving these a formula is obtained for the fractional elongation of the 
cylinder. When the cross-section is bounded by a single closed curve the equation 
satisfied by the potential functions reduces to that obtained by Green and Shield 


for the case of torsion of an isotropic cylinder, and their general method of solution 
applies. 

The problem of the second-order bending of such a cylinder by terminal couples 
is also reduced to the solution of a single boundary-value problem for two complex 
potential functions and the classical boundary-value problem for torsion. A general 
formula is found for the change of length of the line of centroids and the boundary- 
value problem is solved for the case of right circular cylinders. 


1. Introduction 


GENERAL methods of successive approximation for problems in the theory 
of large elastic deformations of isotropic bodies have been proposed by 
Signorini (1), Misicu (2), Rivlin (3), Green and Spratt (4), and Rivlin and 
Topakoglu (5). Rivlin has applied his method to the problem of the torsion 
and extension of cylinders of arbitrary uniform cross-section as far as terms 
of the second order. This second-order problem has also been discussed 
by Sheng (6). 

A solution of the problem of second-order effects in the torsion of isotropic 
incompressible uniformly extended cylinders has also been given by Green 
and Shield (7) in terms of complex potential functions. They have shown 
that complete details of the stress distribution can be obtained once certain 
integral equations are solved. Blackburn and Green (8) have obtained 
corresponding results in their investigation of the torsion and bending of 
compressible isotropic beams. 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 2, 1958) 
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In the present paper the theory developed by Green and Spratt (4) is 
extended to the case of a material which is transversely isotropic with 
respect to a certain constant direction. The results are applied to a homo- 
geneous incompressible cylinder which is transversely isotropic with respect 
to its generators. Complex variable notation is used. It is shown that the 
second-order problem of torsion and extension can be reduced to the 
solution of a single boundary-value problem involving two complex 
potential functions. In the case of a simply connected cross-section this 
boundary-value problem is identical with that obtained by Green and 
Shield (7) for an isotropic incompressible cylinder, and once this is solved 
the complete second-order stresses and displacements are known. Without 
solving the boundary-value problem a formula is obtained for the fractional 
elongation of the cylinder. The problem of the bending of the cylinder by 
couples over the end is also examined as far as terms of the second order. 
The solution again reduces to a boundary-value problem involving two 
complex potential functions together with the boundary-value problem 
for the classical! torsion function. The boundary-value problem is solved 
for the special case of a circular cylinder. A formula is obtained in the 


general case for the angle of twist and fractional elongation of the central 
line. 


2. Notation and formulae 


Theories of finite elastic deformation of aeolotropic materials have been 
developed by Ericksen and Rivlin (9), Green and Wilkes (10), and by 
Adkins (11). The notation of the present paper is based upon that of the 
latter author and of Green and Spratt (4). We summarize the relevant 
formulae here. The points P, of an unstrained elastic body initially at rest 
are defined by a system of general curvilinear coordinates 6*. If r(@', 67, #) 
is the position vector of P, referred to a fixed origin and ¢,, g' are the 
covariant and contravariant base vectors of the curvilinear system 6* at 


P, then - ; ‘ 
. g:=T;, g'=9"8; 9°9n = % 
Jixk = 8i-Bx g* = g'.g*, 9 = Jul 
where g,;, g are the covariant and contravariant metric tensors. Indices 
take the values 1, 2, and 3. A comma and a vertical line will be used to 


denote partial and covariant differentiation respectively with respect to 6°. 
When the body is deformed a point initially at P, is displaced to P and 


(2.1) 


—>> 
we assume for the displacement vector P, P = d an expansion of the type 


(6, 62, 63) = ev(6, 62, 63) + <2w(61, 62, 63)-+..., 
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where v and w are vectors and « is a real parameter which is regarded as 
sufficiently small for subsequent operations to be valid. We write 
d= di, ¥ $75 a” Sn» v= Un s” — v" Sm) W = Un yg = w"S mn: 
(2.2) 
Then the covariant and contravariant metric tensors G;,, G* in the 
deformed body can be expanded in the form 
Fie = Pie e(U4| eA) er(Ws |W Eels OM le) + 
Qik — gi*— e(vi k4 vk |f) —€2(wi k 1 ayk = ele 
ym |kyyi 1k |moyi 
; ; —y™ |kyt| —vk mot) +... 
and the invariants of the deformation 
,=9'G,,, I, = I3g,,G"*, 3 = Gig = |Gy 
can be written as 
f= +- Zev" | +-€2(2w" |, +-v" |*v,|,)+... 
I; ‘+. dev" rt e*(40 |,+ 2v" -v*|,—v’|,v*|, +0" |*v,|,)4 
I, = 1+-2ev"|,+-€?(2w"|,-+- 20" |, v*|,—v"|,v*|,)+.. 
When the material is incompressible G = g and hence we must have 
Zu |, = v"|,v*|, (2.5) 
The stress tensor per unit area of the palieanensl body and the component 
of the body force vector per unit mass, referred to base vectors at the 
point P of the deformed body, are given respectively by 


She tik 9 Wik 
gk — e7'th ¢2p "tk | 


FR = «€F’'+<2 Fk... 
The equilibrium equations are then 
r'tk| +p, F’* = 0, rik) 4 kl rt 4 (ul | rt) | 4 


ig7 


where p, is the density of the undeformed material. 


3. General theory 

We consider a material which is incompressible, homogeneous, and trans- 
versely isotropic with respect to the z-axis of a cartesian coordinate system 
(x,y,z) in the undeformed body. The curvilinear coordinate system @ is 
chosen such that 6* = z and 6! and 6? are functions of x and y. Then for 
a homogeneous transversely isotropic body the strain energy per unit 
volume W is a function of the three strain invariants J,, J,, J, and the 
quantities J,, J; defined by 

I, = 9" Go—1, Ts = (9 Gag —85)(9"* G4 —87). 
For an incompressible material J, = 1 and the formula for the stress given 
by Adkins (11) becomes 
zik — g*O+(1, g*—g''gG,. yy +M*$+- N*L+G*p (3.1) 





SECOND-ORDER EFFECTS 
where 
gH wu 8h Tae ee 
él, al,’ ol,’ 
Mu — MY MS — M27 — M2? — M3 — M3! — NY , (3.2) 
Nu — N12 — N12 — N®— 0, M® = } 
NB — N31 — ]gkG,, NB N3 — jgQ,,, N® — 4G3_} 


Om — 


22W 
al, 


and p is an arbitrary scalar function of the coordinates @'. From (2.3), (2.4), 
(2.5), and (3.2) it follows that 
I, = 3+eJj+..., I, = 3+713+... | 
I,= elgtelt+..., J, = el5+... , 


5 


(3.3) 
where 
== 9)(7"+e), =e 
Ii = us 37 $0 |>vm\35 rs — } (v3 |" -+-v7 |§)(vs|,+-0,/5) } 
®, ‘’, ¢, and # are evidently functions of J,, J,, J,, and J,;, and assuming 


expansions as power series in € to be possible, we obtain from Taylor’s 
theorem 


aw ew 
) = 2w + 2ev3|,{ ——_ wey VE = 204 2en4|,{ ——_ ie 
: bial if Tat). - — erat) 


+- 2evy’|,+-... , (3.5) 
t+ Beye? |g teP{dv(v9| +0 |9)(vg|,+ Urls) +3(0*|g)*+ 
+y(2w*|,-+-v™ |3v,, |3) + 27v,|,(0" |?+-0*|")}+-... 


(3.4) 


\ 


where 


a). <=ee = (OW = (%) 
él, € ° : él, ay J él; ont is el, ia 


2( oc) o2W 5 3 - (Fr 
a\(\@W- 9 y= ——_ . == mney P = amugeennnges 9 
y OT4} =o elt} e=0 OL, 1s) .=o 


(Sz ew 
ng ee le 
él, ol, ' a,el,) <0 


We now substitute from (2.3), (2.4), and (3.5) in (3.1) and note that, 
since there is no stress when the body is undeformed, x is zero and p is 
—(u-+-2c) when « is zero. Hence we obtain 


rik — p'gik 4 (yi |k4 yk |t) 4 dy Mik |, + deem t(yt|k+ yk |?) 
77th = p"gtk + (2703 |,—p' + fvm™)(v8|* + vt |f) + 
+ (ut dam) (wt |* + wk f+ vy, |Rom |t) — 
— (4+ 20)(v8 |! +0! |)(v¥| +0; |") 4+ M*{270, |, (0 |?+-0* |") + 
+8(v3|5)?+-y(2u4|,-+0" |, 0,,/3) + 
+ v(v9\"+ 0" |9)(v3|,+0,|3)} 
L 
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where m@) — m?2 — m2) — m2) — 0, m2 = m@) = m® = m® = 1, 
m®®) — 2, and p’ and p” are scalar functions of the coordinates @‘. The 
repeated indices i, k in terms involving m“ are not summed. 


4. Deformation of a cylinder 

The deformation of an elastic cylinder by forces acting on the ends has 
been considered by Blackburn and Green. We summarize here the relevant 
results. The unstrained body is taken to be a cylinder of constant cross- 
section R, whose generators are parallel to the z-axis where (x,y,z) is a 
cartesian coordinate system. The plane ends are given by z = 0, z = 
We take the curvilinear coordinates to be defined by 


gi == z+ iy, 62 = x—ty, 6 = 2%, (4.1) 


and write (f, f,z) for (6', 62, 68). Then the only non-zero components of the 
metric tensors are 


$2= 4, Gs=1, ge@=2, g@=1. (4.2) 


Hence covariant differentiation reduces to partial differentiation. We note 
that taking the complex conjugate of any tensor is equivalent to inter- 
changing the indices 1 and 2. 

We take the cross-section z = constant of the undeformed cylinder to 
be bounded internally by a number of non-intersecting closed curves and 
externally by another closed curve, the equation of the complete boundary 
being given by G(¢,) = 0 where G is a real function. Then since there is 
no external stress on this surface we have 


7'lk di—7'*k dt — 0, rik di —7"*k dt = 0, {G(l,f) = O}. (4.3) 
The cylinder is assumed to be held such that 
Adi Adi Ad 2 
i a oe Pa at (4.4) 
ez a «al 
at the origin. Then the surface tractions on the end z = / are equipollent 
to a couple $,(eM''+e2M"'+...) and a force g,(eN’'+<?N”'+...) at the 
centroid (f,,,%,,1) of the end z = | of the undeformed body, where 


N= [[73dR, M = —i|{ (fq) dR 


R R (4.5) 
N®=[[7dR, M®=im[{ (f—{,)'"dR 
RK R 


N= [{ {9 +401),,7'"+07|, 71} dR, (4.6) 
iP 
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N= [{ {9 +04|,7°*+0'|,73} dR, (4.7) 
R 

M"! =i [{ {0% vr 4+(Co—L)r"™+0|,7%+0|,7°%)} dR, (4.8) 


R 


Ms = [f {o%r"4(C—La)(r"™+0!|,,7°™-40"|,2'18)} dR, (4.9) 
R 


5. Torsion of a cylinder 

We now consider a cylinder of incompressible homogeneous material 
transversely isotropic with respect to the z-axis, extended and twisted by 
the action of stresses over the ends which are such that N"!, M"!, N”!, and 
M"! are zero. We take « to be such that the force and couple on the 
surface z = / take the values (34+2«+y)eAR and (u+4x)eS respectively, 
where S is the classical torsional rigidity and A is a real constant. Thus 
we have 
N’3 = (3u+2e+yAR, M2 =(p+}n)S, NB=0, M'=O0. (5.1) 
(The choice of these values for VN’ and M’ eliminates almost all ratios of 
elastic constants from the subsequent expressions for the stresses and 
displacements. ) 

We consider the case when the body forces are zero and look for a solution 
as far as terms of O(e?). The solution of this first-order problem is well 
known and we recapitulate the main results here. 

The first-order displacement is such that 


v= —H+ikz, 8 = —-Hl-ikk, =F =eM+HS(+F(O)}, (5.2) 
where f (£) is a function of f. Then the displacement gradients are given by 
vl}, = —}iz, w|/,=0, H|,=—f, &F, =—3f’, lp =A 
M1 Vy\2 = —P—hiz, %4\3= —hid, rg], =4f', slp =A }. (5.3) 
v1 — 0, v/?— —A+2iz, FHP= iY, F=f’, wBF=A 

Introducing these values into (3.6) we see that 
Moa7R— 0, B= (utguy(f'+il), 7° = (3u+2e+y)A, (5.4) 
if we choose p’ to be Au. The boundary condition on the surface G(f,Z) = 0 
ee (F'+itdl—(f’—il)dt = 0, — {G(E,£) = 0}. (5.5) 
We assume that f’ is regular and f+ is single-valued in the open region R 
of a cross-section of the cylinder, and that f’({) and its derivatives are con- 


tinuous onto the boundary except possibly at isolated singular points. The 
torsional rigidity S is defined as 


S = 1+i ff tf’ dR = I- {sv dR, (5.6) 


0, 


- 
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where J is the moment of inertia of the undeformed cross-section about 
the z-axis. 

Combination of (2.5) and (5.3) yields for the second order incompressi- 
bility condition ' 3) 


;= P24 40(¢f'—Z’), (5.7) 


while from (3.6) and (5.3) the second-order stresses are found to be 


a Cw yl “~ Be ra *“¥\9 
7 — (a —2ibf' +04) 20 F+ibh 
3 
78 — aph—p( 252 + 2] 42202 Wt 2olf"—iL)(F' +L) 
1 3 - 
£18 (u}de)( 4 24 Le tAP— BIN) \ (5.8) 
od (27—20—2n-+ Jw) f’+iL) 
7733 p’+(u+kK 24 it) + 


+ (dy + 48+«+-3v+-77—80—5y)A? + 
+(r—p— 20+ Joy f' +i0)(f’ — if) 


The second-order page equations are obtained by substituting from 
(5.3), (5.4), and (5.8) in (2.7). Thus, 





r ” 42. 3 
cp 12 2,,.3 o*u ee 
— (w+ 43x) V2w + («4 zr = 0 
Oa" 22541 2,3 5.9 
op Ou u : . 5.$ 
2 P + pV2w! L die —- 432° =—- + 2if’—f ~ ( 
eC 7 oie 0Cez 


rs 2o(f’f’+if’— iff”- t) = 0 
We look for a solution for equations (5.7) and (5.9) which will make all 
stresses independent of z and 7”! a constant multiple of 7’. Guided by 
(5.1) we assume that 


ws = 1(B—1)A(f+f), (5.10) 
where f is a real constant, and we note that this restriction on 7”!* implies 
that wl = }i(1+26)ALz—}l2? + H(C, 8), (5.11) 


where H is a function of { and Z. On substituting from (5.10) and (5.11) in 
(5.7) and (5.9) and eliminating p” we obtain 


H eH oH ®H 8H - 
4S = 4+ (L’—Z7") J oY (5.12) 
a * e avreE alot? 
By integrating the second of these equations we deduce that 
oH oH f+if , 
: “= (w+ 2o)fif+if+40'(—10'(O}/p, (5.13) 
ral C4 


where ((¢) is a function of { and the expression if-+ if has been introduced 
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for later convenience. On adding to the first equation of (5.12) and inte- 
grating we find that 

c 
H(L,£) = Q2E+4i j tf'(t) dt+o | {f(D}? di(2u)+ 
0 


c 
+ (u-+20)a(2)+ cO’—042i [ fH dt +2itf} [ay (5.14) 
0 


where @() is a function of Z and the final term has been included to simplify 
the expressions for the stresses. 


In order to satisfy the condition that N”' and M"‘ are zero, it is necessary 
to add to (5.10) displacements which, according to the linear theory, give 
rise to stresses on the end z = / equipollent to an arbitrary bending couple 
and force through the centroid of the undeformed plane end. Such displace- 
ments are well known in the classical theory. It is found in the classical 
linear theory that if 
pl = —)K(28—3lz*) + fee? inde + (KOE + Ri g— Kye) + 

+ }(kl2— kL, —kLE,— 2h) 
v3 }(KE+KC—Kl,,—K¢,,)(2?—2lz)+ 
+(h— Akg BRE + BEG t+ bl a2 + HF (0) + F(D}— 
(4+ 3x4 2y)(KEE+ ROL—2K Ug —2K Eb q)/{16(2n-+«)} . (5.15) 
M™ = 0, 

N" = }(34+2«+y){K(I,+1,)+K(1,—1,+ 2i1,,)} 

N’ (3u+2k+y)AR 

M” bi(3u+2e+y){k(L,+1,)+k(L, —I,+21,,)} } 
where J,, J,, and J, are the moments and product of inertia of the un- 
deformed cross-section, A is a real constant, k and K are complex constants, 


+ = im ff Eo — O80) F’+4(8p-+ 2+ y (Ko + Ktlo)— 


R 


then 





-2(4+ 3x + 2y)KEl—(8u+-5«+ 2y) KO} d R/{8(2n+-«)(1,+1,), 
and: F(Z) is the classical flexure function for an incompressible transversely 
isotropic material, defined as the solution of the boundary condition 


{8(2+-«)(F’ +ial)+-4(34+2x+y)(Klg+Kil_)— 
—2(4p4-3x+ 2y) KCl — (84+ 5+ 2y) KO} dg 
= {8(2u+-«)( F’—ial)+4(3p+ 2«+y)(Kl + Kilq)— 


—2(4p4+3x + 2y) KCl —(8. + 5x + 2y) KE} dt 
=0 (G@=0). (5.16) 
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Thus from (5.10), (5.11), (5.14), and (5.15) we have 
w! = —LK(z3—3lz*) + 4(k—{)z?+ 4i(1+ 28)dALz+- 


— +Kl_—KE)(z—l)+iale+ 
—ktl g—klE g@—2hl) + 9A + 
c c 
+ ft oF i to [ {7} dé/(2u)+ 


0 ry 


c 


a(f)+¢0'—2+42i | f(t) 
a 


Cg—Klq)(2*—2lz)+ 
F)+4(B—1) pe f)- 
+ 2y)( (KEE + of — 2K fC ,—2KCEE,,) {16(24+-x)}. 
y substituting in (5.9) from (5. “t and integrating we deduce that 
du{(K0+ Kl—Kl,,—Kl,)( tke +kog—kl —kl+ 2h 
+ 2i( f—f)—Q’ a ah C2} + 4u(if—if+ CZ) 
+o(f'f'+itf’—iff’ —Q’—O'—20)+(20427—p)dA*. (5.18) 


From (5.8), (5.17), and (5.18) we obtain 


— (w+ 20)(Q’ +O’ +222) 
}(3p+2«+y)(KEl,+ Kil ,)— 
Hau 3x + 2y) KE — dy (8u+-5u+ 2y) K+ 
bv-+27—20—2u+Bu+ $Bx)A( f’+i2) 
y){(KG+KE—Kl,—Kl,_)(z—-l)+ 
1 h+kl +k —ke— FQ +(r— _p—o+}y)x 
ur TNS So +4) f—f)— 
— + 3(3u43x+y—4o)lZ 
+(4 ss 1§+«+-3y+ peep ) 





Introducing these expressions into (4.3) and using (5.4), the boundary 
condition on the curved surface becomes 
d{(u+20)(@+60’ + Q)—i(2Qu+x)f(f—f)+(20—n—p)l2Z} = 0 
(G(,Z) = 0), (5.20) 
together with the boundary conditions (5.5) and (5.16) for the classical 
torsion and flexure functions. 
By multiplying (5.20) by { and ({—{,,)?, integrating around the com- 
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plete contour and transforming by Stokes’s theorem we obtain the following 
equations which will be used later, 

[[ (e+ 20)(2'+2')—i(2p+-n)(f—F)} dR 

R 


= (2Qu+«)S—(4o—x)I 
J] —fa){(u t+ 20)(Q' +O" + 200)— 
R 
— (2p +x)[0E+-if—if+ (0—if’)(f—L,)]}} dk = 0 


When the cross-section is simply connected we can integrate the first- 

order boundary condition in the form 
f—f= itt (G(e,0) = 0). (5.22) 

Then (5.20) can be written 

B+00+Q4+0E=0 (G(,2) = 0). (5.23) 
A boundary-value problem of the same form as (5.23) was solved by Green 
and Shield (7) in their investigation of the torsion of an incompressible 
cylinder. They obtained integral equations for Q and w when the cross- 
section is simply connected and found explicit solutions for particular 


cross-sections. These results are immediately available for the present 
problem. 


. (5.21) 


The second-order resultant forces and couples over the end z = / are 
derived by introducing (5.3), (5.4), and (5.19) into (4.6)-(4.9). On using 
(5.5), (5.15), and (5.21) we obtain 

N"! = (3u+2«+y)idlg R+4K(L,4-] )+4K(1,—I,+2il,y)} | 
N”3 = (3u+2«+y)(AR+41)—(2u+o—r+4n—}r)S+ 
+ (48+ 3y+3«+3v+97—60—3yu)A2R 
M"! = (3y-+2e-+y)ALo(l+-BA)R+ WAL, +1,) + 
+K(1y—I,+ ily) +4 [[ Awe S+F VF’ +) — 


R 
~(S—Lal(r—p—o + oy f’ +40) fi) + 
+ 3(3p + 2e+-y)CE+ F(Qn+K)(2f fF’ +ilef’— 
—iff’+lol—C)}} dk 
M"? = (34 +2«+y)ACL, +L) +AS{(u+ dx)B+ 
+ 4v—}hx—20+27—3y} | 
Since these are all zero we can determine h, k, K, and 8. It is found that 
B = {(6n-+-40—47+«—v)S—2(3n4 2x+y)(L,4-L,)}/{(2u+«)S} 
h = {(84-+40—47r-+ 2«—v)S— 2A? R(54+-3y4+6x+ . (5.25) 
+ 6v-+ 187—120—6y)}/{4(3.4-2«+-y)R}—I/(2R) 
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We note that K is zero if either the axis of torsion is the locus of the centroids 
or there is no first-order extension. In this case the expressions (5.15) and 
(5.19) for the stresses and displacements reduce considerably. The angle 
of twist per unit length of the undeformed cylinder then becomes 


e7A(B+-1). 


For an isotropic material, y, 5, x, v, and + are zero and these results 
reduce to a slight generalization of those given by Green and Shield (7). 


6. Bending of a beam by couples 

We next assume that a cylinder of homogeneous incompressible material 
transversely isotropic with respect to the z-axis is bent by couples acting 
on the ends z = 0, /. We choose the origin to be a point of the surface 
z = 0 of the undeformed cylinder. The displacement at this point will be 
restricted to satisfy (4.4). For convenience in the subsequent calculations 


we choose « such that r 


= (83u+2k+y)eR, 
where T is the magnitude of the deforming couple, and look for a solution 
up to terms of O(e?). We begin by giving the main results of the first-order 
problem, which are well known [see, for example, Love (12)]. 
The first-order displacements are assumed to be such that 
vi = $M224+-41(Ml2—Mil,—M&Z,), 
v — —}(Ml+ MC—Ml,—Mlq)z, ne 
where M is a complex constant. Then the displacement gradients are 
y/, = 0, v!|, = Mz 
LiL, —Mfl— it) 
LMCq, v,|35 = 4Mz 
-_ —Mf— Mz), yljl — 0 
- Mz, v3|1 — —Mz 
= 1 Mt, —Mf— a0) ) 





Hence the first-order stresses are found from (3.6) to be given by 
8 an an gO on 7? — }(3u+2«+-y)(MZ,,4+ Ml,—Ml—MNzd), 
(6.3) 


provided that we choose 
= 4u(Ml,4+ME,—Ml—Mb). (6.4) 


We note that these displacements and stresses satisfy the first-order 
incompressibility condition (2.5) and equilibrium equations (2.7) and leave 





SECOND-ORDER EFFECTS 153 


the surface G(f, 2) = 0 free from stress. From (4.5) and (6.3) we see that the 
resultant of the stress on the end is a couple 


(3u+2e+y)eRe® = fi(38u+2e+y)ef{M(L, +1, )+M(L, —I,+2iI,,)} (6.5) 


where @ is the angle between the a-axis and the direction of the couple. 
We deduce from (6.5) and its complex conjugate that 


u — iR{((L+1 eee —I,+2if,,)e- 
pri 2(13 I,) 





(6.6) 
By writing (6.5) as 


[{ (Mf —LoNE—Lq) + M(L—Cq)} dR = —2i Re, 


R 


we obtain the following results which we shall use later : 


| ( (M{Z—ME,,4+MC—Ml,) dR = 2iR(Me-®— Me) 


J } (mete + It, + MM ,— Ml, — IPC —MMAE,) dR eee 
= 2iR(Me-"4- Me) 


From (2.5) and (6.2) we deduce the second-order incompressibility con- 
dition to be 


w' |, MM2+-} (IP? 4+ MM El+ M2?) + 
+ fe(ME q+ Moo) MEq-+iq—2ME—2M10), (6.8) 
whilst from (3.6), (6.2), and (6.4) 


(4 F + Ms) 
C6 


2p” (2 + Mae .s 
C2 
+4(u—20—27)(ME + Ml—Ml,— Mfg)? 


be tee HM 4 Me BME — BMT) Me} 


"+4(2u+ 2b) (252 + Me) + 
+3(y+6+ 2«+6v+ 14r—160—10p) x 
x (Ml+MC—ME,—Mf,)?! 


By substituting from (6.2), (6.3), (6.8), (6.9) in (2.7) the second-order 
equilibrium equations are found to be 





5.” 
cp 


+ (ut Jo V0 + MT) + (0 b( SE+ M2) =0, (6.10) 
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wn” ey ae | - 
Fe + uVrol + de So +2 aa) te 2y+-«)(2M%— Ml ,—MML,)— 
Cc Cz C 


—(20+27+3x+4y)M(MlZ+MC0—ME,,—MEL,,) = 0. (6.11) 


We look for a solution of equations (6.8), (6.10), and (6.11) in which 7” 
(and thus r”*4) is proportional to z and the remaining stresses are indepen- 
dent of z. Since the first-order displacements v!, v® are respectively even 
and odd functions of z we assume that w! and w’ are also respectively even 
and odd functions of z, and thus we write 


wi = Lf, f)z?4-A(f, £), ws = —1MM2+ F(f, Z)z, (6.12) 


where F, H, and L are functions of £ and {. The incompressibility con- 
dition and the conditions imposed above on +r”! imply that 
OL : ¥4 I 59 L — 
4— +1 art? -+MM = 0, 


C4 et 


and hence that 
L = k(k+ipt)- +MMC—M*Z,—MM{,,), (6.13) 


where f and k are real and by constants respectively. The restriction 
imposed on 7”! implies that F is harmonic and can therefore be written as 


F = h+ 4kl,+-3klq—3k0- . Mf? 4 M202 + 
MME g4+ Ml) ME4+-ME)+4B{f(2)4F(}, (6.14) 


where h is a real constant and f is an arbitrary function of ¢, this particular 
form for F being chosen to simplify subsequent results. 
If we introduce these expressions into (6.8) we deduce that 


CH eH = 7 “ 
at t oe ~ Be Te Ma—® BU SAS )+fe(ME4+MC,)?— 


OM ME)+4(3M2224+-2M MCZ+3M22). 
(6.15) 
By substituting for w! and w* in (6.11) and eliminating p” by differentiating 
with respect to ¢ and equating to its complex conjugate, we obtain 
@®H @H 
fn ES = —48. 6.16 
atzae atet? 88 (6.16) 
We find it convenient to integrate this in the form 
BE ag = 2-20) + 4B FSi) + EE 
soe ea (6.17) 
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where Q(Z) is a function of £. By adding equations (6.15) and (6.17) and 
integrating, we have 


{ 
H = &(f)+ 0-04 RO — kel — Bel g— 2ho)— diperl—ap | f(t) dt + 


-§( ME, + Mle) —( ME + Ml) MC+ EME + eM MCE, (6.18) 
where @(£) is a function of Z. Thus we conclude that 
wl = }(2k4- 2160+ M2l,,4+MMl,—M%—MMf)z*+ 
+ (RL? — kel — kel g— 2hl) +H + 60’ —Q— pipet — 


4 


AB | f() dt +4 ML + Mea) — HME + Mog) M e+ 


Ly MMOL 4 EME 
—4M Mz*+-{h—4kl—3k0+ ph0g+3klo+3B(f+F)+ 
+ (MME + MME 4+ Ml,,4- M*{l,— M*t?@— Mt*)}z] 





Since the stresses and displacements must be single-valued and continuous 
we require that Q’+ 0’, 40’+-Bf, and 6+¢0’—OQ—4B fro dt should be 
continuous and single-valued over R. Rigid body deplocmente can be 
added to (6.19) to satisfy (4.4). 
From (6.10) and (6.14) we see that p” is independent of z and hence by 
substituting from (6.19) in (6.11) and integrating, we find that 
p” = pth+ pklq+ pho g— ght — gk —20'— 20 —3( M+ M20) — 
—}(ME,+Ml,)?—4(ME,+ M_)(ME+ Mo)}+ 
+ }(40+-47+ 2x+-y)(Ml+Ml—M~,,—M,)?. (6.20) 
From (6.19), (6.20), and (6.9) we obtain 
rll = p(40' +400" +1 MM +3M2CC— hip?) 
78 = {MM el— 1M*t2— 1M*t?— §(Mlg+Moq)(Ml+Ml)— 
— Bf —Bf—42'—40} + }(2«+y)(ME+ ME—Ml,—Ml,) 
7718 = $B(2u-+«)(f’ +i) 
7°83 = (u+K- by )B(F+F)+ (34+ 2e+-y)ht bebe + del g— tk — 
— hl — 4M? — 4 MC?) 4 Ml + gM t AMM lt 
+4M MEl,4+-4(2y+8+4n+6rv— 126+ 187)(MZ+ ME— 
—ME,,—ME_)*—p{2Q' +20’ +4 M224 $M Mil+ 
+1 M2t?—4( Ml + Mog) Ml+ Ml)+3(Mlo+Mlo)"} 
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The boundary conditions on the surface G = 0 are found from (4.3) and 
(6.21) to be given by 
eG 
‘= 
(Hn -+y)(ME-+ MC —ME,— MC q)*+ pl M MCE — 2 i M¢2)— 


p{4a’ +400" + 4M M04 3.MPCl— hist? 


~3(MZ,.4- MC) MZ+ Ml)—4(0' ++.0')—B(f+f pe =-0 (@=0), 
(6.22 
(f’+il)df—(f’—if)dl = 0 (G=0). (6.23) 


The condition (6.23) is the same as (5.5) and the values of f are known for 
many cross-sections. For later convenience we choose f (0) to be such that 


[J F+P)aR =0 


If we multiply equation (6.22) by Z and ({—Z,,)?, integrate around the 
complete contour, and transform 4 Stokes’s theorem we deduce that 


i ( {16u(Q’+-0’') dR 
R 
( ( {(3u+2«+-y)( MeG?— M*C?— M22l,,-MMK,+ MME, 4+M%X,)- 
‘R 
—4Bul( f’+il)+-10upM M04 2(2«+-y)(MlZ+ Ml—Mf,,— ML,)?- 
— Ol rh +Mf,,)(MZ+Mg)} dR, (6.24) 


[ [ v6u(Z—Z,)(2' 
= 
‘ | | (—f,){uf3(M M64 M2l)\(2—L£,,)+-4MMCl— 32 —3iP22]— 
3u(Ml,+MC,)(3Ml+2MC— ME,)—2Bp[2f+ = (C—f,)( f’ +0f)]+ 
+(2e+y)[2(2— + (C—Lq) MI (—f,) M +(f—Lq) MJ} aR. 
(6.25) 
On equating to zero the real and imaginary parts of (6.24) and using 
(5.6) and (6.7) we obtain 
({ (Q’4+0’)dR 
= 5MM1/8—9(Mf,,4+-M{,,)?R/16+ 
+-i(2«-+-y)(Me-— Me") R/(4p) 
B = (84+2«+y)(Me i@ +-Me®)R (2u8) 7 
We note from (6.5) and (6.26) that 8 is zero when the couple is parallel to 
a principal axis of the undeformed cross-section. 
We can obtain the resultant forces and couples on the end z = / by 


(6.26) 
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substituting in equations (4.6)-(4.9), from (6.1)—(6.3), and (6.21). By making 
use of (6.23) and (6.26) we find that 

N “=0, M%=0 
N"3 = (3u+2«+y)(hAR+4MMI)— . (6.27) 
— }i(12u+ 120—187—6v—2x—y—8) R(Me-— Me) 
M" can be found similarly but for simplicity we give here only the value 
obtained when £,, is zero. In this case 


M"! = hi(3u+2x+y)k(L+1,)+kU,—I,+2i1,,)}+ 
hi [ { {2(8u+2+-y) MC? +-(15p+ 120—187—6yv—8) x 
R a 
x (ME+ Ml) —2Bp(26f+Cf’ + 20f—ill)} dR. 


Since NV"? and VM"! are zero h and k can be determined from these equations. 
The value of h is 


h = i(12u+ 120—187—6v—2x—y—8)(Me-— Me) /{4(3u +-2x+y)}— 
—MMI/(4R). (6.28) 


We note on using (2.3), (6.1), and (6.19) that the length after deformation 
of the locus of the centroids of the undeformed cross-sections is given by 


t t 
VG dz\ " | [1 +ev3 |, +-€2(w|,+ 4v1|, v?|,)+...] dz\ 
0 ‘ 0 


Fe Iz-%6 


1 +-eh+ 4M Ml Fo+48f (La) +48F Ea)]+--}- (6.29) 
We consider as an example the case of a circular beam, given by 
Gl, 2) = @—a*. 
We take 6 = }n, so that we find from (6.6) that 
M = 4/a?. (6.30) 
From symmetry we see that 8 and k are zero, while from (6.28) we have 
h 2(94—4n—2y+120— 18r—6v—8)/{(3p+2«+-y)a7}. (6.31) 
We assume that 
w = M*w, €, Q = M(Q,f+-Q, &), (6.32) 
where w,, Q,, and Q, are real constants. Then equation (6.22) hecomes 
4w, €+- 120, a2f—12Q, —8Q, C+$af— 
—FO4+42Qc+y)(C+OC/p = 0 (@=0). (6.33) 
By writing { as a®/{ and equating separately to zero the coefficients of £3, 
¢, and 1/f in (6.33) to determine w,, Q,, and Q, we infer that 


O = (5p-+4ee+2y)f)/(2a2)— (3. — 2x —y)2/(3a*), 
w= —2(3u+2«e+y)l/at. (6.34) 


Se 
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The corresponding displacements and stresses can then be determined from 
(6.19) and (6.21). 

The magnitude of the applied couple is 

P = (3u4+2«+-y)na*e. (6.35) 

Then from (6.29) and (6.31) we find that the increase per unit length of the 
central line of the beam is 

l’ 2(9u— 4x — 2y+ 120— 18r— 6v—8)T? 

l ra®(3u+2«+y)* 





(6.36) 


The author wishes to thank Professor A. E. Green for suggesting the 
problem and for many helpful remarks, the Department of Scientific and 
Industrial Research for financial assistance, and the referee for his helpful 
criticism with regard to the presentation. 
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SUMMARY 


The mixed boundary value problems with which the paper deals are problems in 
which the plate is clamped along part of the boundary and is either free, or subject 
to specified bending moment and shear, along the remainder. Complex variable 
analysis, and especially the techniques evolved by Muskhelishvili (1) and (2) are 
used throughout.. First the general boundary conditions for the problems are stated 
and reduced to a suitable form and then, using known results for plates with fully 
clamped boundaries, the boundary conditions are reduced to non-homogeneous 
Hilbert problems—the kind of problem treated in Muskhelishvili’s book (2). Special 
attention is given to the case of point loading in the interior of the plate, and to the 
case when part of the boundary is subject to constant bending moments and shears. 


1. Introduction 

Tue theory of the solution of mixed boundary value problems in plane 
elasticity is based on the work of Muskhelishvili and his Russian collabora- 
tors, and results of their work are to be found in the books (1) and (2), where 
further references are given. In these books attention is directed to the 
solution of problems in plane strain and generalized plane stress primarily 
from a general standpoint, although a large number of special problems, 
such as those connected with punches and stamps, are also studied. Green 
and Zerna (3) also consider a few mixed boundary value problems connected 
with punches and cracks. The solutions of two mixed boundary value 
problems in plate theory have been given by Cooper (4). The first of these 
is that of a circular plate clamped along part of its boundary and acted on 
by constant bending moment and shear along the remainder, while the 
second is that of a circular plate with a partly clamped-partly free boundary 
acted on by an isolated normal load at the centre. 

The object of this paper is not to present the theory of mixed boundary 
value problems in plate theory in all its generality, nor to develop an analysis 
which is mathematically rigorous; the mathematical techniques have 
already been rigorously developed in (1) and (2), where they have been 
applied quite generally to plane strain problems in many ways mathemati- 
cally analogous to the problems of plate theory which are treated here. It 
is relevant to remark that the analogy between these two fields, though 
quite strict (for some boundary conditions) is only formal; the knowledge 
of the solution of a problem in one field rarely allows one to write down the 

(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 2, 1958] 
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solution of the corresponding problem in the other (cf. (1), 319 and (3), 253). 
In this paper particular problems of thin plate theory are studied and their 
solutions are either given explicitly, or reduced to a form amenable to 
computation. 

Mixed boundary value problems are considered for plates occupying the 
following regions: 


(i) the upper half-plane ; 
(ii) the interior of a circle ; 

(iii) the infinite region bounded internally by a circle. 

The method may be adapted to apply to a larger class of regions by using 
conformal mapping, but this is not done here. 

The plates are clamped along a part C, of the boundary C and are free, 
or have specified bending moment and shear, on the remainder C, of C. 
There are two main types of mixed boundary value problems, namely, 

(a) C, is free and the plate is subjected to an isolated normal load at an 

interior point ; 

(b) C, is subjected to given bending moment and shear. 


The plan is as follows: after stating the fundamental equations and 
deriving some general results, the boundary conditions for problems (a) and 
(b) are given in a convenient form. Two particular problems of the half- 
plane are then solved and finally the necessary modifications are noted for 
the corresponding problems in cases (ii) and (iii). 


2. Fundamental equations 

Referred to rectangular Cartesian axes O(x,y, Z) with the origin in the 
mid-plane of the plate the transverse displacement w of the mid-plane 
satisfies the biharmonic equation (5) 


Ofw a4 p (2.1) 
é220z2—t—i«S PP’ ¥ 


where 2h is the plate thickness, 2h P = Dis Love’s plate constant, z = x+-iy, 
and p is obtained from the boundary conditions over the faces of the plate, 
which are 


z= 0= yz = 2 over Z = —h, 

zz = O= yz, 22 = —p over Z=h. 
It is well known (see, e.g. (3), 231; (5)) that the general solution of (2.1) is 

w= Wyt+2Q(z) +2Q(z)+ew(z)+a(2), (2.3) 
where Q(z) and w(z) are analytic functions of z in the region occupied by 
the material, and w, is a particular integral of (2.1). 
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3. Boundary conditions 
The boundary conditions on C,, the clamped part of the boundary, are 
w= 0= ewlén ond, (3.1) 

and these may be written (6) 
w= 0= 
Cc 


§ 


infC@ , Ow . 
= fe (+i 7a) on C\, (3.2) 


a 


where @/én and @/és denote differentiation along the normal and arc length 
of C respectively, and the directions of n and s are those of x and y turned 
through the angle a. Using the notation of Stevenson (6), the boundary 
conditions on C, are 


ns = M, nz—énn/és = S ony, (3.3) 
where M = 0 = S when (, is free. 

Equation (2.3) enables the boundary conditions (3.2) to be written in 
terms of the complex potentials; it is desired to transform (3.3) in the same 
way. If F and K are respectively the resultant shearing force acting on an 
are AP of C,, and the complex moment about the origin of the stresses along 
the are, both measured per unit thickness of the plate, then it may be shown 
((3), 232) that 


. P 
F an * | (y, dz — ¥, dz}, (3.4) 
of 
P 
K+izF = 4 {(A+2iF) de—T di}, (3.5) 
A j 


where the notation is that of (5), namely, 
A=ay—yx, T=ayt+yx+2iyy, Y= a+izz. (3.6) 


Equation (3.5) together with the relations 


1 eta, aie = (2A 4 Pe-May Deriay, (3.7) 


dz 
which are valid on the are AP, give 
ref (K+i2F) =e =M on, (3.8) 
Also, from the relations 
nz = 4(Ye-*+ He), n= = (Te-#—Tera), (3.9) 


it follows that 


im 2 (K+icF) = Fie, (3.10) 
ee eee ee 
oF iets TE mm (Yet Tyee) = te, (3.11) 


6092 .42 M 
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and therefore 
ic BD Misesu dee — Onn : ' 
im — —(K+izF) = nz———_- = S on. (3.12) 
és dz es . 
Equations (3.8) and (3.12) hold for both isotropic and anisotropic thin 
plates; they show that specification of bending moment and shear on C, 
enables K +-izF to be determined there apart from an arbitrary combination 
of the form iCz-+-e, where C is a real constant and « is a complex constant. 


It may easily be shown ((3), p. 234, where the notation is slightly different) 
that eee ——— 
K+izF = 2P(1—n){—«Q(z)+20'(z)+w'(z)}, (3.13) 
where x = (3+-»)/(1—7») and nis Poisson’sratio. By using this last equation 
together with (2.3), the boundary conditions (3.2) and (3.3) become 


ZQ(z)+-2Q(z)+(z)+w(z) = 0 on G, (3.14) 
Q(z) +20'(z)+'(z) = 0 on, (3.15) 
KQ(z)—20'(z)—w'(z) = M(z,Z) on C,, (3.16) 
where M (z,Z) satisfies 

—2P(1—7) re : M(z,z) = M, (3.17) 

dz 
—2P(1—1) im — nd M(z,2) = S. (3.18) 

és dz 
In particular, if C, is free (case (a)) then one may write M(z,z) = 0. In 


this case, if 02,(z), we(z) are the complex potentials for the corresponding 
problem when the plate is clamped along the whole of the boundary (the 
displacement function is then the biharmonic Green’s function for the 
region) one may write 


Q(z) = Q,(z)+Q,(z), w(z) = wo(z)+,(z), (3.19) 

where siete _ 
Q(z) +20Q6(z) + (z) = 0 on C, (3.20) 
Q, (2) +2Q4(z)+}(z) = 0 on C,. (3.21) 


The boundary condition on C, may then be written 


KQ, (2) 2.0%; (2) —w (2) = —(«+ 1) Q(z) on Cs. (3.22) 


4. Preliminaries 


Denote the region occupied by the material by S*+, and let S- be the 
complement of S++C. Then, if ®(z) is a sectionally holomorphic function 
(see (2), 33 ff.) with part of C as a line of discontinuity, the boundary values 
of ®(z), as z approaches a point ¢t of C from S+ and S-, are denoted by 
@+(t) and @-(t) respectively. 
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In cases (i) and (iii), when S*+ includes the point at infinity, it is assumed 
that w is at most O(|z|) at infinity. This is consistent with the order of the 
biharmonic Green’s function, which is O(1) at infinity in case (i) (Tiffen (7)) 
and O(|z|) in case (iii) (Dean (8)). Inspection of (2.3) shows that these 
conditions may be written as 

Q(z) = iCz+O0(1), — w(z) = O(1) (4.1) 
for large |z|, where C is real. In case (a), since Qo, wo satisfy such a condition 
it follows that Q,, w, also satisfy this condition. 

In all the problems of type (a) the functions Q,(z), w,(z) are holomorphic 
in S*+ except possibly at infinity (in cases (i) and (iii)). For problems of 
type (6) this last statement applies to Q(z) and w’(z). 

If C and & are real constants and « is a complex constant, then equations 
(2.3) and (3.13) show that the complex potentials 


Q(z) = iCz+e, w(z) = —éz+ik (4.2) 
give w= 0, K+izF = —8P(iCz+e), (4.3) 


and therefore they affect neither the boundary conditions nor the displace- 
ment w. Allsuch terms in the complex potentials may therefore be neglected. 


5. Isolated interior loading of the partially clamped upper half 
plane 
Tiffen (7) has given the complex potentials Q,(z), w(z) for the problem 
of the upper half-plane clamped along the whole of the real axis and acted 
on by an isolated normal load at an interior point. Ifthe load is of magnitude 
R = 32ahP (to exclude a constant from the solution) and acts at the point 
Q (z = 2) then 





O,(2) = —(—%)log(=2), (5.1) 


Z—2y 
z— 
z— 





we(z) = Z9(2—Z) log( 
The function ®(z) defined by 
{ Q(z), zin S* 
\ —2D}(z)—a;(z),  z in S- 


is sectionally holomorphic throughout the whole plane except possibly at 
infinity, where it satisfies 


=) + (2)—Z)z- (5.2) 
9 


O(z) = (5.3) 


®(z) = iCz+0(1), (5.4) 
and transforms the boundary conditions (3.21) and (3.22) into 
@+(t)—@-(t) = 0, on, (5.5) 


«O+(t)+O-(t) = —(x+1)Q,(t) on Cj. (5.6) 
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Here C, is taken to be the union of a finite number of segments of the real 
axis. If the displacement w is to be finite at the end points of C,, which is 
a physical requirement, it is necessary that ®(z) shall be bounded at these 
points. If X(z) is the fundamental solution of the equation 


kX+(t)+X-(t) = 90 on, (5.7) 


which is bounded at the end points of C, and of lowest possible degree there, 
then it is known ((2), 236) that the most general solution of (5.5) and 
(5.6) which is bounded at the end points of C,, is 

(«-+1)X(z) [ Q,(t) dt 


(z) a ~ er X*(pt—2) > P(z)X(z), (5.8) 


Cy 
where P(z) is an arbitrary polynomial. In particular, when C, consists of 
a single segment in the finite part of the plane, the solutions which satisfy 
(5.4) are given by P(z) = iC, where C is a real constant. (If there is a 
solution satisfying the prescribed conditions at the end points of C, and at 
infinity, in the case when C, consists of a number of segments in the finite 
part of the plane, then it is given by P(z) = 0.) In this case, since X(z) 
satisfies (5.7), equation (5.8) may be rewritten as 

; X(z) f Qt)dt oy, 
(z) me Xoera tC*e) 
Cy 
where Cf denotes C;* + Cy (each described so as to leave the corresponding 
region on the left). 
We may take without loss of generality C, to be the segment |t| < a of 
the real axis, for which the function X(z) is ((1), sections 110, 114) 


X(z) = (z—a)*+**(z+a)t-", y = (log «)/2z. (5.10) 


The function h(z), defined by 


h(z) = los( =), (5.11) 


a 
is one valued and holomorphic in the z-plane so cut that no circuit sur- 
rounding 2, or 2, lies wholly in the region of definition; let these cuts be 
made along the lines im(z) > im(z,), im(z) < im(Z)). The function Q,(z), 
which is given by =e 
g ? Q(z) = —(z—2,)A(z), (5.12) 


is then continuous on CJ, and (5.9) may be written 


(z) = eos) ( SOS, axe, (5.13) 


( 


where H(z) = h(z)/X(z) (5.14) 
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. l 
and A= C+ > | H(t) dt. (5.15) 
C3 


The integral in (5.13) is now transformed by considering an integral 
around L, which is the union of the contours shown in Fig. 1. These 
contours consist of a large circle Cp, having radius R and centre the 
origin, Cf, the lines QF, Q’ F’ traversed once in each direction, and small 
circles about the points A, B, Q, and 
Q’. It may be shown that all the inte- 
grals around the small circles vanish as 
the radii tend to zero. Since H(z) is 
holomorphic inside L, satisfies 


H+(t)—H~-(t) = 2mi/X(t) (5.16) 


on QF and Q’F’, and is O(|z|-') for 
large |z|, it follows that the equation 











1 f A(t) dt a 
H(z) = 55 [ Fa (5.17) 
may be written Fic. 1 
1 f A(t) dt dt 
H(z) = 5— — woo 5.18 
«) 2m J t—z | X(t)(t—z) pt 
C ore 


where the second integral is the limit, as R tends to infinity, of the sum of 
the integrals along the lines Q’ F’ and FQ. 
Equations (5.13) and (5.18) now give 


@(z) = (z—z9)h(z)—(z—29)X (z) | x@acgt 4X0 (5.19) 
Fre 


The transformations =o it all an ai 8 (5.20) 
a Il-o a l—wu 
; . dt o® du . 
rield X(z ee eek a -_ 
es te) | xaqaaa | tens Cm Ht oD 
ore r 
A(z) = | Fis (5.22) 
J u—a 
and therefore 
D(z) = (z—2)p(c)+- AX (z) +t Blz—2), (5.23) 


where B is a real constant, [' is the transform of Q’ FQ under (5.20), and 


Ic) = | (1- . (5.24) 
J 


u—o 
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= S 


The terms containing B in (5.23) give complex potentials of the type 
(4.2), and may therefore be neglected. Equation (5.15) may be transformed 
to give 

8 re{A + J,(1)} = 0 (5.25) 
while im(A) may be determined from the condition that there shall be no 
term of the form ik log z, with k real, in the expansion of «,(z) for large |z}. 








Fic. 2 (p 0) Fia. 3 (p 0) 


In fact it may be shown that if the terms in (5.23) are expanded in a Laurent’s 
series for large |z|, then 


(z—29)Jg(o) = al Ip(1) = r eters +0(5) (5.26) 
a z 2? 
‘ 2 
and X(z) = a +ipyt Ps + (5) (5.27) 
where ¢ = 27,(1) +-2(1-+-29/a)Ip(1) (5.28) 
and Pp, = —2(y*+}). (5.29) 
The condition on A then reduces to 


im(Ap,+c,) = 0. (5.30) 


The form of the curve [ may be determined from (5.20); it is given 
parametrically by 


1,8 
Co -- Po- —_ — aX a @ = a, (5.31) 
pot! 
where ef* = —(Z,+a)/(z29+4@), Po = Yre(z,/a) (Fig. 2) (5.32) 


and without loss of generality it is assumed that py > 0. In particular, 
when the load acts at the point z, = iacot(}a) on the axis of symmetry, 
the curve I is an are of the unit circle. (Fig. 3.) 

The function J,(c), which is holomorphic in the o-plane cut along the 
negative real axis (the transform of the free part of the boundary in the 
z-plane) has to be determined by numerical integration. The function J3(c), 
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however, may be expressed in terms of elementary functions as 


: o-1\ 1 
Ig(o) = —[(1— Fi 3) as “|. (5.33) 


where the square brackets indicate that the increment of the function is 
taken along I’. In particular, in the symmetrical case 


I,(o) = 2i{sin ee ee ee (5.34) 


and A=2 [ =e ae. (5.35) 
Equations (5.1)—-(5.3) and sidahieatle (2(z) and w‘(z) to be expressed as 
Q(2) = —(2—29)log{—3) + +(2—2)Ig(o)+AX(2), (5.36) 

w'(2) = - f{log( = =) + flo) ~ 
- 2LIp(o)-+(2 —2)I9(0) + AX"(z) — 2 —- —AX(z), (5.37) 


and these enable the determination of all the stresses in the plate. In 
particular the clamping forces are given by (see (6)) 
yx = 8PreQ'(z), yz = 8Pim Q"(z). (5.38) 

The other interesting quantities are the values of the slope éw/éy and the 
deflexion w along the free edge C,. These may be calculated from (5.3), which 
gives ‘ 
© ww .0w 3 m 
ott * 2{D+(x)—@-(x)} on C,. (5.39) 
Using (5.6), this may be written 

OO x He 41 Q,(2) + (2— 29)I3(0)+AX(z)}, (5.40) 


OL cy 


where the right-hand side is evaluated on the upper side of C,. 


6. Applied stresses on the boundary of the partially clamped upper 
half-plane 

The problem treated here is of type () (section 1). The plate is clamped 
along a part ©, of its boundary and is acted on by prescribed bending 
moment and shear along the remainder C, of the real axis. 

The boundary conditions for the problem are given by (3.15) and 
(3.16), in which the complex potentials satisfy the conditions laid down 
in section 4. The function ®(z) defined by 
{ Q(z), zin S* 


a | —2f'(z)—a'(z), z in S-, 


(6.1) 
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is sectionally holomorphic throughout the whole plane except possibly at 
infinity, where it satisfies 
(z) = iCz+0(1), (6.2) 
and transforms the boundary conditions (3.15) and (3.16) into 
@+(t)—@-(t) = 0 on, (6.3) 
«@+(t)+@-(t) = M(t) on Q,, (6.4) 
where M(t) = M(t,?) is given by (3.17) and (3.18). 

These equations are exactly similar to (5.5) and (5.6) with —(«+ 1)Q(t) 
replaced by M(t); their solution is given by an equation of the type (5.8). As 
in section 5, if C, is taken to be the segment |t| < a, then the solution of 
(6.3) and (6.4) is 

X(z) M(t) dt 


= Drie +1) J Xray t O*). 





In particular, if the bending monet M and shear S are constant, then 
equations (3.17) and (3.18) show that 


—2P(1—n) M(t) = Mt—hiSe. (6.6) 
Again the integral may be determined by contour integration, and the 


condition that w(z) shall include no term of the form tk log z (k real) for large 
iz|, demands that C = 0. Thus 


O(z) = —[(2Mz—iSz*)+{tSz—2(M+ySa)}X (z)]/16P, (6.7) 


and from this equation the stresses and displacements may be determined 
as in (5.38) and (5.39). 


7. Modifications needed for a plate having a partly clamped 

circular boundary 

Consider one such problem, namely that of a normal load acting at an 
interior point of a circular disk whose boundary |z| = a is partly clamped 
and partly free (Cooper (4) has solved the problem for a load at the centre). 
The complex potentials Q,(z), wo(z) corresponding to a load of magnitude 
R = 32nhP at z =z, in the interior of the fully ee disk, are 
(Dawoud (9)) 


bal { 292 1 Zo\\ 
Ole) = —(e—s)llos| ay | -)} 


wo(z) = Z(2 —z4) log! = =e +3( —3\e—), 
0 


L(%p—z)a} 2 
where z = a*/Z,. The function ®(z) defined by 


_ {Q,(z), zin S+ 
= —20)}(a?/z)—a}(a2/z), z in S- 





SOME MIXED BOUNDARY VALUE PROBLEMS 169 


(S* is now the region |z| < a) is sectionally holomorphic throughout the 
whole plane except possibly at infinity, where it satisfies 

D(z) = —20,(0)+O(1), (7.4) 
and transforms the boundary conditions (3.21) and (3.22) into (5.5) and 
(5.6). The general solution of these equations is again given by (5.8). In 
particular, when C, consists of a single arc, which without loss of generality 
may be taken to be |arg(z/a)| < «, then 


X(z) f Q,(t) dt 
O(z) = — rene ft AX(s), 7.5 
) om | Xyt—z) +4 — 
C3 
where X(z) = (z—ae’*)+*7(z—ae-**)t-7, (7.6) 
and A has to be determined by (7.4). 

The function Q,(z) is now defined in the plane cut along part of the line 
joining z, and 2), and the integral in (7.5) is evaluated by a contour integra- 
tion similar to that used in section 5. It may easily be shown that on 
introducing the transformation 

z o—ei® 
$4 ft. 7.7 
a oe'*—] (7.7) 


the function ®(z) may be written 
D(z) = (2—29)p(o)—A _ 3 \(@—%)—AX@), (7.8) 
0 


where J,(c) is given by (5.24) (for a different arc I’) and A is determined from 


Ae-**7(cos «+ 2y sin a)+A = Ipe'2) + IpleF) +201 —e Ie) —[ _ = 


(7.9) 
The curve I’, which is the transform under (7.7) of the above-mentioned 
cuts in the line joining z, and zo, is always symmetrical about the real axis. 
It may easily be shown that if arg z, = 5 and sin(a+5) 4 0 then I is an arc 
of the circle in the o-plane having centre sinScosec(a+5) and radius 
\sin « cosee(a+8)|, while if sin(a+5) = 0 it is a segment of the straight line 
reo = cosa. In particular, if the load acts at the point z) = apy (|p9| < 1) 
on the axis of symmetry, then [ is the are jargo| < « described in the 
anti-clockwise direction, where 


e* = (pp—e*™)/(pp e’*—1). (7.10) 
In this case /2(c) is given by (5.34) with a replaced by ¢, A is real, and 


€ 


Ip(e*) . faentay — — [ Smbr(0—a) 
A” )+4,(e-**) = (aes dé. (7.11) 
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The corresponding problem for the infinite plate bounded internally by 
the circle |z| = a is very similar. The complex potentials Q,(z), wo(z) are 


Q,(2) z—2 — + 
~ 


z 


fale 


Z9 


= Z,(z—2,) log(5— 


“j” «© 
ae (7.13) 
\z Zl} 29 
and the function ®(z) defined by (7.3) (where S+ is now the region |z| > a) 
is sectionally holomorphic throughout the whole plane except possibly at 
infinity, where it satisfies 

@M(z) = iCz+ O(1), (7.14) 

and transforms the boundary conditions (3.21) and (3.22) into (5.5) and 

(5.6). Using the same notation as for the circular disk, the solution is 
given by (7.5), in which A = iC and 

X(z) = (z—ae'*)'-*7(z—ae-**)i +47, 


A contour integration again allows ®(z) to be written as 


(z) = —(2—29){I3(¢ \—log? =| 4-AX(s ), (7.16) 


where J,(c) is given by (5.24) with 8 replaced by £, and in which I is an 
are similar to that mentioned previously in this section. The constant A 
is given by equations similar to (5.25)-(5.30); in fact 


re|4 —Ip(e-**) + } log 20\ == 0, (7.17) 


“0 


im) A(1—e*7(co8 a5 2y sin «))+ Jp(e"*) —Ip(e-**) + =0(1 eta) Ig(e's) = 0. 
(7.18 


In particular, if the load acts at the point z) = a@/py (|p, < 1) on the axis 
of symmetry then I is the same arc of the unit circle as in the case of the 
circular disk, and - 

€ 

* ¢X9-a)_ eos }(0—a) 


asi tee 7.1! 
2 sin }(@—a) 40 + 108 po. ae, 


. 


—€ 
The problems (6) for plates with circular boundaries may be solved in 
a similar way, but this will not be attempted here: one such problem has 
already been treated by Cooper (4) (see section 1). 





SOME MIXED BOUNDARY VALUE PROBLEMS 


8. Conclusions 

It is found that the numerical computations involved in finding the 
clamping stresses (5.38), and the deflexion and slope on the free edge 
(5.39), are not unduly heavy. In fact, in the most interesting case, when 
the load acts on the axis of symmetry, it is found that the three integrals 
I,(c), Ig(o), and I,(c) given in (5.24), (7.8), and (7.16) are all closely linked 
together, and one integral of each type may be obtained from one numerical 
integration: some numerical results have been given by Cooper (4). 

It is sometimes claimed that for each solution of a problem in elasticity 
theory by complex variable methods there is a corresponding solution using 
only real analysis. This may be true for plates having fully clamped 
boundaries, although even there the complex analysis is remarkably lucid. 
But for the mixed boundary value problems treated in this paper the theory 
of sectionally holomorphic functions is essential. This method, which may 
also be used for plates having fully clamped boundaries, has two obvious 
advantages: first it dispenses with the tentative approach, which is often 
no more than judicious guesswork, and secondly it immediately furnishes 
a proof of the uniqueness of the solution under certain specified conditions 
(see Gladwell (10)). 
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SUMMARY 

Previous parts of this series have been concerned with the solution of Laplace’s 
equation in three dimensions. In this paper the methods already developed are 
extended in order to include in their scope the solution of stress problems in three 
dimensions; from the relaxational point of view this requires the simultaneous 
satisfaction of four Poisson-type equations in four dependent variables. The paper 
firstly derives these four governing equations, with corresponding boundary con- 
ditions which express given values of the boundary tractions, and then demonstrates 
their solution for a particular example of stress distribution in a cube. 


1. In the general three-dimensional stress system there are six components 
of stress at any point; these stress components can be expressed in more 
than one way in terms of three stress functions so that the three equations 
of stress equilibrium are automatically satisfied. Two such ways of ex- 
pressing the stresses are well known, and are due respectively to Maxwell (1) 
and Morera (2). We have found Maxwell’s expressions much more suitable 
for the purpose of deriving governing equations capable of solution by 
relaxation methods; when the stress functions are denoted by x, ¢, and y, 
and with orthogonal coordinates x, y, and z, the six stress components are 


Oy od ~ é*y 
oy? T ez? lear Oydz 
ey am ary 
é2z2 * ax?” ~ @zbat 
a My - a 


ae a= — 
oat Gye? = Y= Srey 


rm 0 





We shali find it necessary also to introduce a fourth unknown, the principal 
stress sum 0, given by © = +9742, (2) 


which is well known to satisfy the relation V?0 = 0. 
(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 2, 1958) 
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2. The equations governing the four functions, x, ¢, y, and ©, are derived 
from the equations of strain compatibility; these number six in all, three 
of the type: ne ne . 

. é € 
Sw | 2 (3) 
dz? © dy*® 
and three of the type: 
o Mee o Rye, ee, Cay) (4) 
Cyoz = Ox ex ey oz 
in a customary notation (3) for the strains. When these strains are ex- 
pressed in terms of the stresses, and hence by (1) and (2) in terms of the 
stress functions, equation (3) becomes 


e e 
j{O—(1 +0) V¥9} + —(O—(1+0)¥%9} = 0, (5) 
cy~ C2 


and equation (4) becomes 
- {Q@—(1-+0)V2x} = 0 
eyez “” 
where o is Poisson’s ratio and where 
weet 
dz? © Oy* © dz* 
3. There are three equations typified by (6) and these can be immediately 
integrated to yield 
O—(1+oa)V?x = falz,x)—fe(2, 9) 
O—(1+s)V*d = fox, y)—fiy,) |, (7) 
O—(1+o)V*%p = frly,2)—Sal2, 2) 
where the functions f are each arbitrary functions of two of the coordinates. 
Substitution into (5) shows that 


Of, _ fe 
tat — dy?” 
from which it follows that a function X(y,z) exists such that 


eX ex 
A = ey?’ A == ea’ 


It follows similarly from the other two equations of the type (5) that 
functions ®(z,x) and Y(z, y) exist such that 
ap er er 
f= ==> fs = 53° and fe= Far (9) 


Cz 


(8) 
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By virtue of (8) and (9) equations (7) may be rewritten 


a al 
@—(140)V8y = 0_& 


éa? a? 
wr Fi 
oy? dy? 
eX #&©O 


ez® saz? 


QO—(1+e)V24 = (10) 
and O—(1+o)V¥%s = 


We now introduce functions a(y, z), 8(z, x), and y(z, y) which are solutions 
of the equations 


(l+o)V2a = X 
(l+o)V78 = © }, (11) 
and (lto)V4y = 


and thereby transform (10) into 





6*B o 
@ = (1+0)V%y4+ 56 — 
te) (x ' Ox? a4 
. on 

@ = (140)¥1(g427_ 22) |. 12 

(140) (54+ mt (12) 
A2 a2 
and o= (+o)¥y4+5—S) 
Cz Cz 


4. Now it can be shown by substituting into the expressions (1) for the 
stress components that the additions respectively of (0?8/éx*—é*y/éx*) 
to x, of (@*y/ey*—é*a/éy”) to ¢, and of (@%a/éz*—e?8/é2) to #, make no 
change to those stress components. In (12) we can therefore regard 
(6*8/ea* —0*y/éx*) as being included in y and soon. The equations (12) there- 
by simplify to become V¥y = @/(1+0) 


V*4 = O/(1-+0) }. (13) 
and V4 = ©/(1+0) 


It can easily be verified that also 
v20 = 0, (14) 


and the four equations (13) and (14) comprise the governing equations in 
their simplest form for the functions x, ¢, %, and ©. We are not aware that 
the equations governing a three-dimensional stress system have previously 
been solved as a result of being expressed in this way; they are not only 
of marked simplicity but they are eminently suitable for solution by 
relaxation methods. 
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5. We do not contemplate in this paper a discussion of the satisfactiont 
of boundary conditions on a general boundary but instead propose to show 
how any conditions (of specified tractions) may be applied at the plane 
boundaries of a rectangular parallelepiped. Axes of coordinates x, y, and z 
are taken parallel to the edges of the parallelepiped, as shown in Fig. 1. 


























H G 
E F 
z 
paidiite a a y 
RR EES Fite 
S----- D ae <p ah OS fe ie ETE va Cc 
A B 
Fic. 1 


The conditions which hold on the face BCGF are first investigated: 
specification of the surface tractions implies that the stresses zz, ry, and xx 
are all known over this face, in general as functions of y and z; knowledge 
of the first of these gives ard 


aa, Ss ae, 
Oz0X 
which may be integrated to provide 
ed PP 
oe 2D ans 15 
A | Fede +asly) (15) 
Knowledge of the second shear stress zy yields similarly 
~ he 
ieee J ty dy +1,(2). (16) 


(The lower limits (z = 0, y = 0 respectively) in the integrals in (15) and 
(16) are chosen only for convenience in making those integrals definite.) 
The manner in which the functionst g,(y) and r,(z) are determined will be 
explained later in section 7; for the present we note only that (15) and 
(16) are normal gradient boundary conditions over the plane face BCG F 
for the functions ¢ and @ respectively. 


+ By substitution from the expressions (1), conditions specifying given tractions can, 
of course, always be stated in terms of the stress functions x, ¢, and 4, whatever be the 
form of the boundary surface. 

{ Ten similar functions occur in the corresponding conditions on the other faces of the 
cube; they are displayed in the statements of those conditions in (19), (20), and (21). 
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For the third stress function y equation (2) is used as the governing 
equation on the face BCGF; on substituting for yy and zz from the ex- 
pressions (1), this equation becomes 

ex ex — — od Oy 


Sy 5? wx — = Dmn2? 
dy? © a2? 0u* =x 


(17) 


this is a two-dimensional Poisson’s equation which serves to determine y 
over the face BCGF. There remains to state a boundary condition for the 
function ©; we shall explain, in section 6, how the first of (1), viz.: 


: = - (18) 


fe = T+ 8, 
ey? ez? 


can be employed for that purpose. 
A set of equations similar to (15), (16), (17), and (18) holds over the 
face ADHE: 





ad , a 
ap | dz — Uy) 
i) 

r y 

ous Re oe 

== xy dy —r,(z) 

: (19) 
< e* * 
44d eee 
oy” a alan 

and salen a aed 
dy* © 02? } 


By the same methods four equations can be derived over the faces 
ABFE and DCGH: 


— ( yzdz+p,(x) on DCGH | 
ex 6 
cy ae 


-- [ # dz —p,(x) on ABFE 
0 





— [ zy dx -+r,(z) on DCGH 











Opp é (20) 
oy : 
— | éydx —r,(2) on ABFE 
Fj 
a2 a2 a2 a2 
M4 EE = Of 
éz2 ' da cy? ey? 
ay am 


and yy = 


oz? * da? 








ing 
eX- 


8) 


he 
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Likewise four equations over the planes ABCD and EF GH can be obtained: 





f z 
—f Zxdx +q,(y) on EFGH | 
ob _) 
az x 
an | fx da —q,(y) on ABCD 
0 
y 
—f yz dy +p,(z) on EFGH 
Ox _ ri (21) 
az y 
-|# zdy —p,(x) on ABCD 
0 





Gx® * dys — Gat at 
~ aha @ 
and z =o ay ) 


6. We next explain how, in principle, these boundary conditions are 
employed to determine simultaneously the functions x, ¢, ¥, and 0. A 
relaxation lattice of the kind described in Part I of this series (4) is drawn 
throughout the parallelepiped with lattice-lines parallel to the coordinate 
axes, and values of © are firstly guessed at the points of this lattice which 
lie on the surface. Corresponding values of © at all internal points are then 
found by relaxation subject to the governing equation (14); this is a straight- 
forward boundary-value problem. Next, values of the stress functions x, 
¢, and & are guessed at all points and then modified simultaneously by 
relaxation, each subject to its appropriate governing equation and boundary 
conditions; thus y is determined according to the first of (13) as governing 
equation and to (17), the third of (19), the first of (20), and the second of 
(21), as boundary conditions on the various faces of the parallelepiped. 

When values of x, ¢, and % have been determined consistent with those 
previously determined of 0, the latter are modified by adjustment of their 
boundary values and the cycle of computation repeated. The adjustment 
of @-boundary values is made in order to satisfy the conditions (18) and 
the last of (19), (20), and (21), on the varions faces. Thus on BCGF the 
values of zx at the points in that face are calculated from (18) by substi- 
tution of the values of % and ¢ already determined. We know, however, 
what zx ought to be at these points (this stress being one of the given 
components of boundary traction), and we therefore adjust the value of © 
at each point in BCGF, increasing it by the deficiency in our calculated 
value of zx there. 
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7. Special reference must be made to the satisfaction of the boundary 
conditions along the edges of the parallelepiped. At a point 7’ (Fig. 1) on 
a typical edge BF’, six conditions must be satisfied, viz. (15), (16), and (18), 
together with the first, second, and fourth of (20); whilst the three governing 
equations for x, ¢, and on BF have been stated in (17), the third of (20) and 
the third of (13). Of the boundary conditions, the first of (20) is used in 
conjunction with the y-governing equation (17); the condition (15) is used 
in association with the ¢-governing equation, the third of (20); and the 
conditions (16) and the second of (20) are used with the %-governing 
equation, the third of (13). There remains to discuss the two conditions, 
(18) and the fourth of (20), only one of which is required, as was explained 
in section 6, to adjust the value of © at the point 7. Two equivalent 
possibilities are open to us; we may satisfy (either) one of these remaining 
conditions by adjusting ©, and then satisfy the second by choosing an 
appropriate value for either r,(z) or 7,(z), as the case may be. In order to 
illustrate the procedure, suppose we choose to use (18) in order to adjust 
the value of @; then to satisfy the fourth of (20) we must give an appropriate 
value to r,(z) in (16). A similar choice at the points U’, V, and S (Fig. 1) is 
now determined for us: at U we must use the fourth of (20) to adjust © and 
hence (18) to adjust the value of r,(z); in turn the fourth of (19) must be 
used to adjust © at V, the fourth of (20) to modify the value there of r,(z). 
Finally, at S, the fourth of (20) is employed to adjust ©, the fourth of 
(19) to modify r,(z). Thus consideration of the conditions along the four 
edges parallel to the z-axis, together, enables all of the four quantities r to 
be adjusted. By similar means, consideration of the corresponding con- 
ditions along the other two sets of parallel edges enables all the quantities 
p and q to be modified. 

Minor modifications of the foregoing technique take place at the corners 
of the parallelepiped, where some simplification results from the fact that 
© is known at starting. 

It is not thought to be necessary in this paper to quote finite-difference ap- 
proximations to the various governing equations and boundary conditions. 
Such replacement in finite-difference terms and the subsequent definition of 
residuals, elimination of fictitious values, construction of relaxation patterns, 
etc.,are now standard procedures in the application of the relaxation method. 

8. At first sight it might seem that the determination of a solution 
according to the methods outlined in sections 6 and 7 would be a lengthy 
and arduous process. In fact it is not so, for the reason that the solution, 
to a particular problem, for the functions x, ¢, and % is not unique. (For 
we have already seen in section 4 that y, ¢, and % may be taken to include 
the functions (€°8/éx*—é*y/éx*), etc., respectively; and the functions a, £, 
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and y are themselves merely required to be some solutions of the equations 
(11), hence each is therefore arbitrary to the extent of an additive plane- 
harmonic function—a of y and z, B of z and 2, y of x and y.) 

The existence of the consequent family of solutions for x, 4, and ¢ (all 
yielding of course identical values of the stress components) lends itself 
ideally to the relaxation method. The computation is much reduced by 
leaving x, 4, and & with a (clearly considerable) degree of arbitrariness. 
Restrictions could be imposed in order to fix these stress functions uniquely: 
by doing so we might render the numerical process more precise and its 
explanation easier to describe; but it is a greater advantage to forego the 
fixity, and so to ease the burden of the calculations, whose aim becomes 
merely to discover any one solution of the family. It is in the nature of 
the relaxation method that the fewer the constraints which are applied the 
better, that the less we are compelled (or choose) to specify about our unknowns 
the easier is it to determine a sufficient solution. 


9. Asanumerical example we treat the case ofa uniform cube (sufficiently 
represented by Fig. 1) whose surface is unloaded apart from two equal and 
opposite compressive forces P applied at the centres, J and K, of two 
parallel faces. The solution was computed firstly on a coarse relaxation 
lattice of side L/4, L being the length of an edge of the cube. By virtue 
of the eightfold symmetry, results can be sufficiently presented on diagrams 
showing only one-eighth of the whole cube, and Figs. 2 give the solutiont 
obtained on a second, finer, lattice of side L/8: Fig. 2 (a) records values of 
the stress-functions x and ¢, presented as multiples of 0-0001P; Fig. 2 (6) 
gives values of the principal stress sum ©, as a multiple of 0-016 P/ L*; whilst 
Figs. 2 (c) and 2 (d) respectively show values of zx and yy and of #2 and Zz, 
also as multiples of 0-016P/L?. 

In the computation leading to this solution it should be recorded that 
a priori restrictions, or constraints, were imposed which made the functions 
x, ¢, and & unique. The restrictions actually chosen were that x should be 
zero on EH and ZA, that ¢ should be zero on EA and EF, and that % 
should be zero on EF and FH. It was not initially realized that these 
constraints, which analytically appear to simplify the solution, in fact 
render its numerical determination more arduous. The calculations pro- 
ceeded without undue difficulty, but we nevertheless feel impelled again 
to assert that when a solution permits of the inclusion of an arbitrary part, 
an advantage is thereby conferred which should not be sacrificed in the 
unhelpful interests of orthodox simplification. 

+ Values of ¥, 2, and Zy are not included in Figs. 2 since these may be deduced, by 


virtue of the problem’s symmetry, from those of ¢, 7y, and 2, respectively, by reflection 
in the plane EFJK (Fig. 1). The value zero was taken for o in the calculations. 
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In conclusion we feel that it is desirable to draw attention to the fact 
that, whilst the methods of this paper should prove satisfactory for three- 
dimensional stress analysis throughout regions enclosed by orthogonal 
plane boundaries subject to given boundary tractions, on the other hand 
these methods have not yet been developed to deal with the more general 
case of a curved surface. That development is clearly desirable but it must 
await, as a preliminary, a successful method of satisfying a normal gradient 
boundary condition in three dimensions. 

Grateful acknowledgement is made to the D.S.I.R. for a grant which 
enabled this work to be undertaken. 
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SUMMARY 


The method of Comrie and Bower for subtabulating a function by use of bridging 
differences is generalized to cover all sub-intervals where the procedure gives exact 
pivotal values. Two other methods are given; one uses leading differences and the 
other corrections to all odd orders of differences. The merits of the three methods 
for use with an automatic computer are discussed. It is shown that any desired 
accuracy can be obtained even when the pivotal values are not built up exactly. 
The method can be applied to all subtabulation ratios including those that are non- 
integral or even irrational. Various worked examples are given. 


1. Introduction 

THE problem discussed in this note is that of subtabulation by finite differ- 
ence methods of a function for an interval of the argument other than that 
for which it is already tabulated. Comrie (1) and Bower (2) have described 
how a function tabulated at a given interval can be tabulated at a smaller 
interval by a method essentially involving only summation. However, the 
method described by Comrie and by Bower is only applicable when the 
new interval is an integral sub-multiple of the old and furthermore when 
that interval is simply related to the radix of numerical notation. (For 
example, with the usual decimal notation, the interval ratio has to be of 
the form 2-*5~° where a and b are positive integers.) Making use of the 
fact that an automatic computer can be programmed to work to any desired 
accuracy and so the required number of guarding figures can always be 
carried, the method has been extended to subtabulation ratios which are 
not related to the radix of notation, or integral, or in the final extension, 
even rational. 


2. Synopsis 

The leading differences and two types of bridging differences used in fifth 
order interpolation for an integral number of sub-intervals are tabulated 
in section 3. The relative merits of the methods using leading and bridging 
differences are also discussed in this section. In section 4 expressions are 
given for the accumulation of rounding errors in those cases where the 
process is not an exact one and so the pivots are not built up exactly. 
(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 2, 1958] 
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Section 5 considers the case when the number of sub-intervals in the 
original interval is non-integral. Section 6 gives an example of some com- 
plexity; a function is given at an interval of 1° in the argument and is 
subtabulated at an interval of 0-001 radians. 


3. The leading and bridging differences for fifth order interpolation 
for an integral number of sub-intervals 

We follow Comrie in denoting the central difference operator for the 
original interval by A and for the new interval by 5. If Az is the original 
interval at which the function F(x) is tabulated, the new interval is 6Azx 
(0 < @ < 1) and we require the values of F(z,+-s0Ax) where s is a positive 
integer. It is convenient to denote F(x,+-s0Ax) by F,, 9. 

Using Everett’s interpolation formula with sixth and eighth differences 
thrown back to the fourth, we obtain the following table for leading 
differences. 

TABLE 1 





Aves 2 2 MS (AM). 


Bsa | 8 f,0(1 —6)(2—8) | 450(1—6*)(4— 8?) 
540 ib 9( 1 — 0)(1+ 78) —}0°(1—6*) 
r+ie — }0°(1—36) }6°(1 — 56%) 
8S a9 ch 26° 

5° 40 & 























In this table, M! = A!—0-207A8+-0-045A8, 
(where A* < 15,000 and A® > 200), (3), and 
(AM*),,, = Mt,,—Mé. 


If A® > 15,000 but A® < 300,000 and A? < 30,000, simultaneous throw- 
back from sixth and eighth differences to second and fourth may be used (3). 
In this case Table | is modified by writing M? for A? and (AM?),,, for A? 


r+ 
where M? = A?—0-01312A8+-0-0043A8, 
M4 = A4—0-27827A$+-0-0685A8, 
and (AM?*),,, = M?,,—M?. 


r+1 


The five sixth order bridging differences (see (1) and (2)) are given in 
Table 2. (A description of the method is also to be found in Hartree’s 
Numerical Analysis (4), which may be more readily accessible.) 

Here (A?M*), = (AM*),,,—(AM*),_,. 

If the simultaneous throwback defined above is used A‘ in this table and 

in Table 3 becomes (A?.M?),, where 


(A?M?), = (AM?),,,—(AM?*),_,. 
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The sixth difference may be kept zero everywhere and when the pivotal 
value has been reached corrections can be made to all odd orders of differ- 
ences, The expressions for these corrections are given in Table 3. 


TABLE 2 





Mt— As (A?M*), 
5o 29 = S899 40(1 — 6*) 73541 — 6*)(4— 6) 
Big = Big | —§0(24+-6*) | — pA 1 — 0*)(8 + 130%) 

ae 01 — 6) Ho 9(4 + 508+ 1104) 














TABLE 3 





Correction to| M*—A$ (A?M*), 


bp +t 40(1—6*) |450(1 —0)(4—6?) 
ay ~— — $6°(1—6*) 
8 95 


r+ 














It is always possible to correct only the alternate differences at the same 
level but it must be remembered that in this case the true differences of 
the subtabulated function are never formed as they are when sixth order 
bridging differences are used. 


3.1. Choice of method for bridging from one interval to the next 

These methods are applicable to many types of machines but high-speed 
computers are particularly useful owing to the large quantities of data 
involved and the simple repetitive nature of the procedure. 

The choice of method depends on the value of 6 and the computer to be 
used. If @ is of the form 2-¢5~°, where a and 6 are positive integers, for a 
decimal machine, or 2-¢ for a binary machine, and the time for a multi- 
plication or division is large compared with the time for an addition, 
subtraction, or shift, it is probably quickest to use sixth order bridging 
differences. On the other hand, the organization of the programme is 
easiest when leading differences are used and this advantage may outweigh 
the disadvantage of the greater time taken to evaluate them. 

In all methods the highest order difference involved is (A?M*), which, 
being of order ten, requires five function values before the point at which 
subtabulation commences, No initial conditions are inserted ; the procedure 
starts with the first function value given and, taking all earlier differences 
to be zero, builds up to the second. The subtabulated values of the first 
five intervals are meaningless. If 6 is related to the radix of the machine 
as above, and if sufficient figures are retained, the pivotal values are built 
up exactly and this is a good check on the accuracy of the calculation. 
The unit used should be such that & is an integer. When the subtabulated 
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values are produced they must be rounded off to the accuracy of the original 
data. 

If @ is not related to the radix of the machine, the choice of method 
depends on the accuracy and not only on the speed. If @ is rational it is 
possible to avoid the build up of rounding-off errors by multiplying the 
original values and all formulae by 1/65 and dividing the subtabulated 
values thus found. This method is slower than the unmodified leading 
difference method as it requires one more multiplication and 1/8@—6 more 
divisions per interval of the original table without improving the accuracy. 
The accuracy of the methods of section 3 is discussed in the following 
section. 

TABLE 4 
A A* yr - i 


575 

1465 1952 

4345 5869 

13094 14633 


36476 


4. Subtabulation when 1/@ is an integer but not simply related to 
the radix of the machine 


An example of subtabulation to sevenths (6 = }) is included to illustrate 
the difficulties that arise in this case. 


4.1. Numerical example 


We require F(x) for x = 5(})7 for the function given in Table 4. Leading 
differences are used for the interval (5,6) and then sixth order bridging 
differences are used to obtain the values for the interval (6,7). Since @ = }, 


we have 1/65 = 75 = 16,807 < 100,000 and so we carry five extra figures. 
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Using M} = A*—0-207A$+-0-045A8 and differencing, we obtain: 


x Mét (AM, +t (A?M), 
5 57 
16 
6 73 
170 
7 243 
872 
8 III5 


Using Tables 1 and 4 we obtain the leading differences 


55,, = 1110549, 88, = 114310, 88. = O-18112, —3¢, = 0-02564, 
88,, = 0-00095. 


Using Table 2 we obtain the bridging differences: 
58, = 58, = —0-05498, 
58, = 58, = 0-20616, 
8$ = —0-29321. 
Working with the fifth decimal place as unit, we have 


TABLE 5 
F he 5 & 
22500000 
1110549 
23610549 114310 
1224859 
24835408 132422 
1357281 
26192689 153098 
1510379 
27703068 176433 
1686812 
29389880 202522 
1889334 
31279214 231460 
2120794 
6 33400008 
6 corrected 33400000 2378638 
6} 35778638 296887 
2675525 
38454163 334481 
3010006 
41464169 377134 
3387140 
44851309 425856 
3812996 
48664305 481657 
4294653 
65 52958958 545547 
4840200 
7 57799158 
7 corrected 57800000 
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The rounded-off subtabulated value must agree with the given pivotal 


value; the exact value is used for continuing the process. 


4.2. Analysis of error 
If we denote the error in 5¢ by ef, the error in the subtabulated value 
F.,,9 due to errors in the leading differences is 


.. a. 
e(F.49) = = (“et sn (1) 


5 oat | 


, — k\ . ' ig or . 
for k < 1/@ and where | is the binomial coefficient which is defined as 
J 


zero for k < j. When bridging differences of the type given in Table 2 are 
used, the error in the pivotal value F,,, is 


S (Ue 


— » 


s=—2 


S 


, 6 58 eee 
and 5.4 = 8, Gad ™ Gut 

When bridging differences of the type given in Table 3 are used, the 
error in the pivotal value F 


r+ 


» is given by putting k = 2/6 in equation 
(1) increased by a second term, viz. 


~ (1/6 , ; . ‘ 
S | (correction to jth order difference at first pivot) (3) 
—~ ) 


j=1 


and reduced by the correction made at the first pivot. Similarly, for the 
third and subsequent pivots. 

Applying the formulae of equations (1) and (2) to the example in section 
4.1 the predicted error in the second pivotal value is 84118733, which when 
rounded off gives 842, the actual error. 

In general, if » is the maximum error in any difference, p the radix of 
the machine, and ¢ the number of guarding figures, 


and 


The magnitude of the maximum error in F,,;,9 due to leading differences is 


5 
E . /k ao eel , 
€(F..49) <7 2 (") < 35 (H -5k4*+- 2543+ 5k? + 94k) 
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If s = (p/0)+q and F,, 9 is built up solely from the sixth order bridging 
differences given in Table 2 starting at 58_,, and 1/6 = N, we have 


Ls —Nj+4 —Nj 
«(F...9) = Zz ((° a )+( 5 Net s-s0t 


j=0 
_Nj _Nj41 _Nj42 
+((" TA) +( as JJeteso+(" Jet (6) 


The maximum of ¢ is of magnitude } or }— 1/4 of the last figure according 
to whether 1/@ is even or odd. Let this be y’. Then 


(Maximum error in F,, | 
Pp a ™ y. 
; ‘ s—Nj+5 Ae s—N)j 
J< Ze HD 


j=0 


<¥> ((Meta+9)_ (Mia) 


j=0 


j=0 k=s—Nj 


12? 
< > Wi+al(Nj+4)+15(Nj-+9)*+8}. 
~ j=0 


We consider the error at a pivotal value (when g = 0) since the maximum 
error increases with j and q. 
\Maximum error in F,, ,9| 


15N% 


' 2 
[oF (+ 1t2pt+2p—1) +4 3 


7’ 


24 


<= 


8N 
Pp+1t+Fp(p+))| 


'N ’ . 
< J P(p+1)[N4p(p+1)(2p?+- 2p— 1) +-45N*p(p + 1) +48]. 

It will be seen that for a long table of values, a large number of guarding 
figures is required. For such values of 0, i.e. those not simply related to 
the radix of the machine, it is preferable to use the leading difference 
method in which the number of guarding figures is completely determined 
by 6. It should be noted that however small the error due to rounding off 
the differences, the error in the rounded off subtabulated value may still 
be 1 in the least significant place. 

5. Subtabulation when 1/@ is rational but not integral 

We are given F_,, F_,, %, Fi, F,.--, and we require F,,, F,,.9, Fy,+426s-++» 
where a, < 6. Let k@+¢ = 1 where k is a positive integer and ¢ < 6. If 
%y < >, we have ay +k <1, 
but ap t+(k+1)6 > 1. 

Let ay t(k+1)0—1 = oy, 
then «, is the first point in the interval (x,,2,) for which a subtabulated 
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value is required. Similarly, if «, is the first point of the form a)+-p@ (p an 
integer) in the interval (z,,z,,,), then we have 
Oe = a y+(k+1)0—1 = a, ,+0—¢ (if a. < 4); (7) 
or x, = aytkO—-l=a_,—d (if «_, > ¢). (8) 
Suppose a, is the first value of «, for r= 0, 1, 2,..., s—1, 8,... which is 
greater than ¢. Then a, = ee blh~ds 4, 


a 5 
and so s> —. (9) 
Instead we take U4, = a&,—h = a +80—(8+ 1). (10) 
Every interval contains either k+-1 or k points at which the value of the 
function is required. Clearly 


r+a, 


F.... = F+ a, Ang —F (I —a,)A?— (1 —at)A3, ,+4 
+55 (1—aF)(2—a,) M$ + = (1—a?)(4—af)(AM4),,, (11) 
and the leading differences at F,,,, are found as follows. We here denote 


the forward difference operator at interval @ by D and use E for (1+-D). 
Then 
no 


- 


Sn ’ nen 7 — ob, ] , _ OF 


— b, Mrs ... b, n+8 € 
a, > ("os a a OF hin +5)6- (12) 


s=0 
Thus these leading differences can be obtained very simply from those 
given by Table 1. Expressions for these can be given in terms of @ and a, 
but it is simpler to evaluate the terms in Table 1 and then use the above 
formula for 87... s sng. 

In this case, as when 1/6 is not simply related to the radix of the machine, 
the round off-error will accumulate much faster using bridging differences 
than using leading differences. To check the calculation it is necessary to 
take the values of F.,, .,9 for p = k—9,..., k—1, k, and extrapolate, using 
Everett’s formula with modified fourth differences, to obtain F,,,. This, 
rounded-off to the accuracy of the original data, must agree with the given 
function value. 


6. Subtabulation when 1/6 is irrational 


In this case we need to consider the number of figures to take for 6. If 
we take @ and «, to the same number of decimal places p, then the error in 
either is less than 4x 10-”. Considering the errors in 


2 3 4 5 
Snrtees Sp scpthO> OF 0+ 6» 5F + 40> 5F a +20 and BF + up + 40 
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due to the errors in @ and «,, and summing these to obtain the error in 
F...,.x9 from these causes, we find that the term in A,,, dominates the 
expression. This term is less than $x 10-?(1+-k)A,,, which must be less 
than 4. This gives a simple means of determining p for any particular @ 
and F. 


6.1. Numerical example 

We are given F(x) for z = 65°(1°)75° as in Table 6 and we require F(z) 
for x = 1-170(0-001)1-186 radians. For simplicity, we use modified second 
differences M? in Everett’s formula; this we may do as A* < 1000. These 
are given by M? = A?—0-184A4. 


TABLE 6 
F(x) A A* 
2144507 
101530 
2246037 
109815 
2355852 
119235 
2475087 
130002 
2605089 
142388 
2747477 
156734 
2904211 
173473 
3077684 
193169 
3 3270853 
; 216561 
3487414 
244637 
3732051 


The expressions] for the leading differences are obtained from equation 
(12), which becomes 


oP, Xp + 3nd = 88, not ott hn eno ae (a —O)88t hn a0 
6 & 


The value of 5,,,, is obtained directly from Everett’s formula with modified 


second differences. Now 6 = Tone = 0-05729 57795 13082 to fifteen deci- 


is 


r+a, 


mal places and since k is the integral part of 1/0, we have k= 17. In 
Table 6, 18A,,, < 107, so we let p = 7. Then (5) 

6 = 0-0572958 
and 7 = 1-17 radians—67° = 0-0360620°. 


5092 .42 oO 
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The results of subtabulating the function at an interval of 0-0572958 of 
the original interval are given in Table 7. 


TABLE 7 
x F 
radians ; rounded off 
(67°) 2355852 2355852 
4128-7842408 
2359980°7842408 2359981 
65851126816 
2366565°8969224 31°2054801 2366566 
6616° 3181617 
"2150841 31°4567699 
6647°7749316 
2379829°99001 57 31°7080597 
6679°4829913 
"4730070 31°9593495 
6711°4423408 
"2106393 
6743°6529801 
2399964 5683279 32°4619291 2399965 
6776°1149092 
24067 40°6832371 32°7132189 2406741 
6808°8281281 
2413549°5113652 32°9645087 2413550 
6841°7926368 
2420391°3040020 33°215798 2420391 
6875°0084353 
2427266 3124373 33°4670883 2427266 
6908°4755236 
2434174°7879609 33°7183781 2434175 
6942°1939017 0°2512898 
2441116°9818626 33°9696679 2441117 
6976°1635696 0°2512898 
3°1454322 34°2209577 2448093 
7010" 3845273 0°2512898 
3°5299595 34°4722475 2455104 
7044°8567748 0°2512898 
24621483867 343 34°7235373 2462148 
7079°5803121 
"186 2469227°9670464 2469228 


To check the subtabulation it is necessary to extrapolate from the values 
for x = 1-181, 1-182,..., 1-185. Everett’s formula with modified second 
differences is used, taking r = 1-183 and @’ given by 

68 = 67-+a4,+160+06’. 
This gives 6’ = 3-823890 


(since for convenience of interpolation we need 6’ in terms of #, not measured 
in degrees). 

This gives F(68°) = 2475087-0294545 which when rounded off checks 
with the given value. 

The function in Table 6 is in fact F(x) = 10* tan zx taken from (6). Com- 
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paring the rounded off subtabulated values with the values obtained from 
the tables (5) taken to seven figures, it is found that there is a discrepancy 
of one in the values for z = 1-174 and 2 = 1-184 radians. It is impossible 
to avoid such errors when the same number of figures is retained as in the 
original data. To obtain one more figure in the subtabulated values, it is 
necessary to take at least one extra figure in the original function values. 
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THEORY OF THE SPINNING BALLOON 
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SUMMARY 
Differential equations of motion are obtained for a flexible uniform string rotating 
with a constant large angular velocity about a fixed axis when the string has a con- 
stant velocity along its length and air drag per unit length is proportional to a fixed 
power of the normal component of air velocity and acts in the direction of this 
component. For steady state curves the equations are solved, a formal solution being 
obtained for the case of zero air-drag, and an expansion for finite air-drag. 


1. Introduction 

THERE are many dynamical systems which can be represented to a good 
approximation by the above differential equations and this is especially so 
in textile processing e.g. in ring spinning, cap spinning, unwinding, rewind- 
ing, ‘doubling’, and so on (and practically all yarns go through several of 
such processes). The angular velocities are of the order of 10,000 revolutions 
per minute and thus only the envelope of the whirling yarn is seen. This has 
the appearance of a balloon whence the name spinning balloon is derived. 
However, the actual curve adopted by the yarn does not lie in one plane and, 
in this respect, the system is more complicated than that of centrifuge 
spinning (the theory of which has been given earlier (1)). Moreover, in 
centrifuge spinning, air-drag is essential for stability whereas in the 
spinning balloon it is usually a minor factor, while the great rotational 
velocities show that gravity can safely be neglected. 

The properties of the spinning balloon are of considerable importance 
and several solutions have been put forward based on various and, usually, 
appreciable simplifications of the problem (a review of these has recently 
appeared (2)). However, the equations of motion of the present paper are 
more accurate and more general than any of these, while the solutions 
obtained are exact for zero air-drag (involving only an elliptic integral 
of the first kind), and, for finite air-drag, the expansions can, by taking 
sufficient terms, be made as accurate as desired and would seem to cover 
most practical cases. A simple numerical method of solving the equations is 
described which gives sufficient accuracy for most purposes and the 
asymptotic solution (which most cases appear to approach) is also given. 

From the differential equations of motion and their solutions the main 
properties of spinning balloons can be deduced. These are discussed in the 
final section of the paper and illustrated by curves computed from the 


solutions. 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 2, 1958] 
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2. Notation 
The z-axis is taken as axis of revolution, distance along the curve is 
denoted by s, taking ds/dz to be positive, while differentiation with respect 
to s is denoted by a dot. The velocity of the string along its own path (the 
‘winding’ velocity) is denoted by v and taken positive when in the direction 
of s increasing. The angular velocity is denoted by w radians per second 
and is positive when right-handed about the z-axis. 7’ denotes the tension 
in the yarn at a typical point P on the curve having cylindrical coordinates 
(r,6,z) with respect to axes rotating with the curve. The mass per unit 
length of the string is denoted by m and the air-drag per unit length taken 
to be k times the component of air-velocity normal to the curve (the other 
component being tangential) raised to the power ¢ and to act in the same 
direction. The values of k and { can be obtained from the curve of air-drag 
coefficient against Reynolds number calculated for a cylinder of the same 
diameter as the yarn (3). 
Now the tangent to the curve at P is the vector 
(7, 6, 2) (1) 
in the direction of the radius, the perpendicular to the radius and the 
z-axis, and the z-axis itself respectively. The air velocity at P relative to 
the instantaneous curve is, thus, 
(0, —rw, 0). (2) 
If the direction of this air velocity makes an angle y with the tangent then 
it follows from (1) and (2) that 
cony ad, (3) 
Hence the normal component of air velocity has a magnitude V,, where 
V 


n 


= rwsiny = rw(1 —r262)h, (4) 
It can also be seen that the direction of this normal component (which lies 
in the same plane as the air velocity and tangent and is perpendicular to the 
latter) is (Fr/{1—r26234, —{1 129214, ér6/{1—r26?}4). (5) 
3. Equations of motion 

We first equate the moment, about the axis, of forces on an element 5s 


at P to its rate of change of angular momentum. The combined moment 
of the tensional force and the air-drag (see (4) and (5)) is 


{d( T'r26)/ds—kV§ r(1 —r26?)#} 8s. (6) 
The angular momentum of the element is 


m 8s{r2w+-vr?6]. 
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Hence, since the winding velocity v is constant, the rate of change of this 
momentum is 


m 88v di r2w+ vr?6|/ds. (8) 
Thus we obtain the relation 
d{(T’—m v®)r29}/ds—kVE r(1— r262)4 = mv d(r*w)/ds. (9) 
Integrating, and writing ds = dr/7, we find, using (4), that 
(T— mv*)r26 = mvwr?+ kus | {rb+1(1 —7292)C +0 /7) dr. (10) 
Now, considering forces and acceleration in the z-direction, we have 
d(T'2)/ds+-kV§ 2r8/(1—r262)t = mv d(v2)/ds. (11) 
Integrating, we find that 
(T—mv*)s = K—kwt [ {r5(1 1262) er /#} dr (12) 


where K isaconstant. Again, considering forces and accelerations along the 
tangent we find that T'4+-mestvd = 0. (13) 
The most convenient form of the integral of (13) is 


T—mv* = 4mw*(h?—r?) (14) 

where h/ is a constant. 

We now generalize coordinates and write 

R= rjh, Z =z2/h, S = s/h, (15) 
pe = 2v/(wh), A = 2kh(wh)?-, x = 2K /(mw?h?). (16) 
However, since dR/dS = #, we shall continue to use the symbol ¢ and, for 
similar reasons, ré and z. The basic equations of motion (14), (10), and (12) 
then become 


T —mv? = lmw*h(1— R?), (17) 
(1— R®)rd pR+-n/R, n=A [ f RE+N( | — 292) KC +) rk dR, (18) 
(1—R?)i=y—p, p=a | { R&(1—r262) Dr 9/7) dR. (19) 

Now, since p24 792452 — (20) 

it follows that 

(1 — R2)2F2 — (1 — R2)2—(pR+n/R)?—(x—p)?. (21) 

Thus, by integrating any two of the equations (18), (19), and (21), we can 

obtain a complete solution. 


The following rearrangement of (21) gives an alternative form which has 
certain useful properties. Now, from (18) we have 


2un-+-n®/ R242 [ (n®?/R%) dR = 2 | {(n-+-p.R®)/ R®} dn 


= 2) | ((1 — R2)r6 R&(1 —r262)KC +7} dR. (22) 
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Again, from (19), we have 
p?—2yp = 2 i) (p—yx) dp = —2A | {(1— R22 RO(1 — 262) 8 /P)} AR. 
(23) 


Substituting for 2* from (20) in (23), and combining with (22) we find that 
(21) becomes 


(1— Re)? = (1— R*)*—p? R*— x22 [ (n2/ RS) dR—g (24) 
where q = 2X | *RO(1—r262)i-2( 1 — R2)r AR. (25) 


The problem can be reduced to that of finding ¢ and ré from (18) and (24), 
for z does not appear explicitly in these two equations. 


4. Zero air-drag solution 

We now solve the equations of motion when k, and hence A, are zero. 
[t follows from (18) that ” is a constant N, say, and p a constant which can 
be taken as zero without loss of generality. Instead of the parameter x 
it is more convenient to take one of the values of R, R, say, for which # is 
zero. It then follows from (18), (19), and (21) that 


x* = (1— Ri)?—(wR, + N/R)? (26) 


and we can write 
R*(1— R*)*7? = R*(1— R*)?—(pR?+ N)?— x2 R? 
= (R}—R*){A+ B(R}— R*)—(R2—R*)} (27) 
where A = R}{2(1— R?)+p?}— N/R}, B= 3R?—2—n. (28) 
Hence, from (19) and (27), we have 
s|* = dZ/dR = Ril RY— RA + BURR R4)—(RI— R24}. (29) 
If Z = Z, when R = R,, and if we write 


Ri— Rt = Ag, (30) 
then (27) becomes 
R*(1— R*)*?7 = A*g?(1+ Bg?— Ag‘). (31) 
It follows from (30) that 
Agdg = —RdR; (32) 
hence, integrating (29) with respect to R, we have 
Z = Z,— x\g—Bg?/6+- (4A — B*)g*/40-+-...]. (33) 


It follows from (27) that if A is positive (negative), R, must be a maximum 
(minimum), and, hence, from (30) we see that g is always real, while (31) 
shows that we should take g to have the same sign as fF. 

We can in fact integrate (29) exactly, for writing 


R?— R*? = bcos*%é, A+bB—b? = 0 (34) 
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we find that (29) becomes 

Z = Z,¥ x{bl(A+b*)}! [ dé/[1—{b?/(A +-b%)}sin2e}. (35) 
However, it is more convenient to express (34) in terms of the minimum 
value of R, R, say, which occurs when = 0. Hence 


Ri—(R,)? = 6. (36) 
It follows from (36) and (28) that 
R?( R)*{ R2+-(R,)?—2—p?}+ N? = 0, (37) 
and (35) becomes * 
Z = 2, xCt | dé/[1—C{R}—(Ryy}sin*€}, (38) 
where C = 1/{2+p?+ R?—(R))*}. (39) 


Thus R, can be found from (37), and hence C from (39) and, finally, (38) 
can be evaluated. Alternatively, ifthe minimum R; is taken as parameter, 
R, can be found from (37) and then C from (39). Again, by virtue of (32), 
it will be seen that, if in (26), (27), (28), and (29), R, is replaced throughout 
by Rj, (26), (27), and (29) are still true (N.B. A and B, but not y, will have 
different values, however). 

The case N = 0 is of considerable interest, for (37) then shows that Rj 
is zero and such curves therefore intersect the axis. In this case, too, we 
have from (18) and (19), 


d0/dZ = h0/ = r6/(R2) = pR|(Rx) = p/x. (40) 
Hence 6 = (u/x)Z+ constant. (41) 


Non-zero values of N can arise even when the air-drag vanishes. Then 
annular ‘restraining’ rings concentric with the axis, on which the yarn 
rubs, are used, and, for the balloon curves below such rings, N will have 
non-zero values. In such cases # can only be expressed in terms of € as an 
elliptic integral of the third kind. However, @ is not of much importance 
for, though end conditions have to be considered in any given problem, 
this involves the tangent and hence 6, but not, in general, @ explicitly. 

The above shows that there are three essential parameters ., NV, and either 
R,, or R,, which, when specified, determine the shape of the balloon curve, 
while the parameter h determines size. 


5. Finite air-drag solution 

Here we have two additional parameters A, (see (16)), and ¢. We now 
solve the equations of motion (18), (19), and (24) by expanding ¢, r6, in 
powers of g (see (30)). Thus we write 


R21 — R22 = A%g{1 +a, g+a.9?+...], (42) 
n= N—ng—n.g—.... p= — x Pigt+peg*+.-.] (43) 
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where the «;, n,;, and p, are to be determined. In fact n, and p, can be found 
immediately from the integral expressions for n and p in (18) and (19) and, 
from (24), we can find a, in terms of n, and p, (if we use (21) instead of (24) 
very involved expressions arise containing n, and p,, though when all can- 
cellations are made we ultimately obtain the same results). Equation (24) 
also justifies the form we have adopted for 7? in (42) though this will become 
apparent later. The expressions for the higher «,, »;, p; are somewhat 
lengthy as is usual with non-linear equations but a fairly economical pro 

cedure is now given. From (42), (43), (30), and (21) we see that 


R*(1 — R?)2(1 —r26?) 
= (1+ P19 + Pag? +...) Rj —Ag?)+ A*g*(1 +04 9+...) 
= (xR, )(1+8,9+89?+-...), (44) 
where @, == 2p,, 8 = pi+2p,—A/R?+ A?/(,xR,)?, 
83 = 2(p, Po+ps)—2Apj/ Ri +a, A?/(xR,)’, 
84 = P+ 2p, py+2p,—(A/Rj)(pj+ 2p_) + A*a9/(xR,)?.... « (45) 
Thus we have 
(1— R®)-1RE-1(1 — 292)kC-Y — (yR,)Sf 1 +8, 9+... 4-9/1 — RF+4+- Ag? 
= (xR) 1/1 — RD) Y(1 +49 +t g?+...), (46) 


where 

t, = $(C—1)s,, t. = $(C—1){8,+(C—3)s}/4}—CA/(1— R}), 
tz = 3(C—1)[83+ 3(C—3)s,{82+ (C —5)s?/12}]—4$2(¢—1)s, A/(1— Rj),... . (47) 
Further, it is convenient to write 


(1+ag+agg?+...)-# = 14+8,9+B9?+..., 


where 
B, = —a,/2, B, = 3a2/8—a,/2, B,; = —5a3/16+-3a, a,/4—«,/2,..., 


(48) 
and also to write 


(l+-t9-+t.9?+...)(1+B, 9+...) = 1+ujg+u,g*+..., (49) 
whence u, = t,+8,, U, = t,+f,t,+8s,... . 
It then follows from (18), (40), (49), (42), and (32) that 
m = M(xR,)+4/(1— RP} f (14+-8,94+...(1+ug+...(—dg) (50) 
and, hence, we have 
nm, = A(x R,)S+1/(1— R3)}, n/n, = 8,4+-4,, 3ng/n, = 8.44, 8, +Up,... . 
(51) 
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To evaluate p it is convenient to write 


(l+u,g+...)(uR?+N—pAg?—n, g—n,g?—...) 
= (wRI+N\1+0,9+rg2+...), (52) 
whence 
(uRi+N)(u,—v,) = m, (uRi+N)(ug—v2) = ny U,+n.+pA, 
(up. R{+N)(us—v3) = 2, Uy+(HA +M)U, +3... « 
Then we have from (19), (46), and (52) 
pix = AxRy)(wRE+N)/(1— YY [ (1+e.9+--.+pig+--M—dg), 
(53) 
whence ' 
Py = A(x Ry)P-wRF+N)/(1— RP, 2p2/P, = +71, 
3Ps/Py = Va +P, %y+P2-.- . (54) 
It can be seen from (51) and (54) that m, and p, are immediately determined. 
We shall now express a, in terms of n, and p, using (24). Now it follows 
from (26), (28), and (30) that 
R*(1— R®)?—p?R*—2yN R®— y2R?+2R? [ (N2/R%) dR 
Ry 
= A*g?*(1-+ Bg?—Ag*). 
Hence, from (24) and (42), we have P 
A%g*(1+a,g+...) 
A*g*(1-+. Bg?— Ag*)-+ 2 Ft [(n?—N2)/R3] aR-q), (55) 
R, 
From (30), (32), and (43), we have 


R3 J ((n2—N?2)/R°} dR 


> Ri ine (n19-+M_9?+...)?—2N(n,g+...)}/(Ri—Ag?)*|(—Ag dg) 
> Am, g>+m.g*+...], (56) 
where 3A Rim, = 2Nn,, 4A R?m, = 2Nn.—nij, 
5A Rim, = 2Nn,—2n,n.+2Nn, A/ R}.... . (57) 
Again, from (25), (46), (49), and (52) it follows that 
q = 2X { [ R2(1 — R2)%¥2 RE-1(1 — 262) Rr /{ R21 — R2)r}] dR 


2A(xR,)P Mu Ri+N)A 


R31 — RE 
x | g(1+a,g+...)(1+0,9+...)(1—Ag?/ Ri)“ —dg) 
— (A*/ R7)(q, 9° +-92g'+-...), (58) 
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where 3q, = 2A(xR,)-"pRE+N)/(1— RPS, 
442/39, = %+%, 5qs/(3q,) = a+, 0,-+0,+A/Rj,.... (59) 

Hence, from (55), (56), and (58), we have 

4 = 2M +G, % = B+2m,+q2, a = 2m3+G3—(A/Ri)(2m, +), 

ag = 2my+-qy—A—(A/Ri)(2m_+ G2), «5 = 2m,+q,—(A/ Ri)(2m,+495),-.. - 

(60) 
When a, is known, then ¢,, w,, ,, and m,, p, can be determined ; hence, so 


Can a, then fy, Ug, v2 and so on. 
It follows from (19), (42), (48), (32), and (52), that 


Z= | (2) dR = x f (lt+p9+...)1+B,9+--.M—dg) 


= constant — x[9+(p,+8;)9?/2+...]. (61) 

All the above expansions have limited radii of convergence and, for 
example, will not converge if g > g’ where #(g') = 0. However, we can 
extend the solution in the following manner. Assuming that the turning 
point of R next to R, is R, (the zero air-drag value of (37)) we can calculate, 
to a first approximation, n and p by using the expansions (43) first about 
R, and then about R|. Thus we obtain as first values for n and p, at R,, the 
values N, and P, say. To obtain a more accurate value of R, (remembering 
that * = 0 at this point) we solve the equation 

0 = RX1—R2)*—(wR*+N,)*— RXyx—P,)?. (62) 
From the appropriate root Rj, say, of this equation, we can expand again 
to obtain better values of N and P and then repeat the process until 
sufficient accuracy has been obtained. This whole process is repeated at the 
next turning point of R, and so on, and a complete solution can thus be 
obtained. 

We can, if necessary, find a better approximation to the next turning point 
than R by solving the equation 

L+ayg+a.g?+a,g*+a,g* = 0, (63) 
for rapid methods are now available for solving quartics (4). 

In most cases the first two or three terms of the expansions for m and p 
give sufficient accuracy since air drag is usually small compared with 
centrifugal force. In any case the main effect of air drag is to deflect the 
curve until it tends to lie tangentially to the local air velocity, when the 
air-drag normal component will naturally be small; in fact by analogy 
with calculations made in the special casey. = 0,{ = 2-0, (5), balloon curves 
should ultimately approach the asymptotic solution, 


ré — 1, f‘=-i= 0, 3R = (w?+3)*—p. (64) 
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Here the yarn lies in a circle, the tension exactly counterbalancing centri- 
fugal force, and the normal component of air velocity is zero. 

In most practical cases a numerical solution of sufficient accuracy can be 
simply caleulated by taking the zero air-drag values of * and ré and, using 
them in the integrands of the integrals of (18) and (19), obtaining a good 
approximation to the correct values of rf and z and, hence, by (20), #. If 
in this process, we replace dR/f by dZ/z and 2d R/r by dZ we can obtain 
second approximations by repeating, treating the first (zero air-drag) 
values of Z as a parameter and using the better values of 7, r@, and 2. 

Direct numerical integration of (18) and (19) by modified Adams— 
Bashforth methods can be used, but the divergence of the expansions in g 
when this is greater than a certain value seems to be reflected in large values 
for higher differences at certain points. Again, further difficulties arise at 
the turning points of R for, here, since ¢ is small, loss of significant figures 
occurs upon evaluating * = {1— r292 2?}t (while sometimes the ‘predicted’ 
values of r@ and 2 are such that (r6)?+-2? is greater than unity), and some 
trial and error is needed to obtain consistent values. The curve A = 2-0 of 
the figure was computed this way. It should be noted, however, that, in 
practice, though A can be as high as 2-50, the corresponding value of R, is 
usually less than 0-25. Our case A = 2-0, R, = 0-5 was in fact chosen for 
illustrative purposes. 


6. Discussion of some balloon properties 

The main features of spinning balloon curves are of great importance 
and we discuss, in this section, some qualitative and semi-quantitative 
properties which may be deduced from our solutions. These are illustrated 
in our figure by curves of R against Z and of X against Z when A = 0, 1-0, 
and 2-0. 

To present the curves more realistically we also give the reflection of the 
R, Z curve. Thus the R, Z curves give the envelope of the rotating balloon, 
and the X, Z curve the instantaneous curve as seen by an observer at a 
distance when a steady illumination is provided plus a stroboscopic flash 
operating once per revolution at the instant when the tangent to the 
instantaneous curve at the apex lies in the plane of the paper. In the 
case A = 0, the X, Z and R, Z curves coincide. 

Most spinning balloons start or finish on the axis and the properties of 
such curves are the most important. In the absence of air-drag such 
balloon curves are somewhat sinusoidal in appearance especially if R, is 
small (in fact as R, + 0 the spinning balloon approaches the case of a 
circularly vibrating string). When R, is not small the curves are less 
sinusoidal (see the case A = 0 in the figure). The distance between successive 
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zero values of & is, for small values of R,, approximately m/v2 = 2-22... . 
This can be seen by letting R, > 0 in (26), (28), (34), and (35). When R, is 
not small this distance is less and, in the case A = 0, R, = 0-5 of the figure, 
the value is just under 1-85. 

Small values of air-drag do not appreciably alter the R, Z curve from the 
zero case except at the points where, in the latter, R = 0. Here R is not 
zero but a finite minimum (which is proportional to A to a good approxima- 


. a ER 



























































A=t0 A=20 


Fic. 1. Examples of spinning balloon curves when { = 2-0. The R, Z curves are 
continuous, the X, Z curves broken 


tion, provided other parameters are the same). In the case of curves starting 
on the axis the second such minimum is approximately twice the first, 
the third is three times tke first and so on, until the minimum becomes 
comparable with (he maxima. In the illustrated case, A = 1-0, at the first 
minimum R = 0-157, at the second R = 0-302, while the corresponding 
values of Z are 1-82 and 3-51. At the first maximum R = 0-500, Z = 0-93, 
and at the second maximum R = 0-518, Z = 3-74. The distance between 
successive minima (or maxima) is slightly smaller—though not greatly 
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different from the zero air-drag case (as can be seen from these figures. 
However, more pronounced air-drag (e.g. the case A = 2-0 of the figure) 
makes these effects more marked. At the first minimum of the case A = 2-0, 
R = 0-291 (just under twice that for A = 1-0) and Z = 1-68. At the next 
minimum, however, F is only 0-472 and Z = 2-98. Subsequent maxima 
and minima are closer together and both approach the asymptotic value 
of R (which, in this case of » = 0, is, by (64), 1/V¥3 = 0-577...). Simul- 
taneously ré and z approach | and zero respectively, and this tendency can 
be observed from the increasing inclination of the X, Z curve to the Z-axis 
and the decreasing values of 7. 

When air-drag and winding velocity are zero (i.e. when A = p = 0) but 
N is finite, the spinning balloon curve has a smaller distance between 
successive maxima (or between minima, which are now finite) for the same 
value of #, than the zero air-drag case, though it is not easy to demonstrate 
this simply. This effect can be qualitatively observed from the curves 
A = 1-0, 2-0 as Z and, consequently, n increases. 

The effect of winding velocity may be qualitatively different however. 
When positive it can be seen from (18) that it has a similar effect to that of 
increasing n. But, in the important case of unwinding, v (and hence ,) is 
negative and it tends to reduce the effect of air-drag. It can in fact cause 
changes of sign in 0. Usually in unwinding the yarn is pulled off a stationary 
cylindrical package through an eyelet on the cylinder axis above the 
package. At the balloon apex (the eyelet itself) since R + 0, by (18), n > 0. 
Therefore for small but non-zero R, ré = pR+O0(R+1). Hence 70 is nega- 
tive here. At the point of leaving the package it must, however, be positive. 
So since n increases as we go down the curve from the apex it must increase 
sufficiently (i.e. the balloon curve must be sufficiently long or of sufficient 
radius or both to make @ positive at the package). Ginzburg (6) has in fact 
observed that 6 initially increases and then decreases but has not given any 
explanation of the phenomenon. Provided that the values of k and ¢ are 
chosen to give the right air-drag at the balloon maximum radius (R,) no 


great inaccuracy will occur if the wrong value of { is chosen. Thus, suppose 
k,, ¢, and k,, ¢, are chosen so that 































k, Vi = ka Ve (65) 

where JV, is the normal air velocity at R,. Then approximately V, = wR, 
a we have : 

and we have k(wR,)o = ky(wR,)S. (66) 


Now the main effect of air-drag is embodied in the term n,g (see (51)), 


*} TA y 
where ny = MxR,)F+1/(1— RYE. (67) 
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But A = 2kh/(wh)?-$ and y = (1—R}), if n and p are reasonably small 
(see (16) and (26)). Thus we have, as our two estimates of n,, 


2k, ho lwo -2(1 — RR) REH, 2k, h&-lw'-2(1 — R2) RE, (68) 


By (66) these are the same. Some numerical tests have indicated that 
minima might be underestimated by as much as 8 per cent. if £ = 2-0 were 
used instead of f = 1-5. 
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NOTE ON A SOURCE IN A ROTATING FLUID 
By J. WALTON (Aeronautics Department, Imperial College) 
[Received 14 February 1957] 


SUMMARY 
The effect of introducing an irrotational source into a rotating fluid has been 
studied by Barua. In this note, results are obtained which are based on a different 
hypothesis from that of Barua’s. The relation between the maximum flow rate 
which the source achieves and the speed of a long wave on the surface dividing two 
régimes of flow is also pointed out. A numerical solution is found for the shape of 
this dividing surface by means of relaxation methods. 


1. Introduction 
THE flow due to a point source on the axis of a fluid possessing solid body 
rotation has been investigated by Barua (1). The fluid is taken to be 


ROTATING 


Y 





AXIS OF 


SYMMETRY 


inviscid and the flow from the source is irrotational in the sense that there 
are no vorticity components and that the azimuthal velocity is zero. 
The discharge forces its way in both directions along the axis of rotation 
where the pressure is lowest; in the final steady state the irrotational flow 
is confined to a cylinder of radius r, except for a ‘bulge’ of radius r, near the 
source (cf. Fig. 1). Barua assumes that circular elements originally within 
this cylinder expand and conserve their circulation as the flow develops. 
This assumption, however, leads one to suppose that rotating fluid initially 
lying between the cylinder r, and the contour of the final bulge (Fig. 1) must 
be swept away to infinity while circular elements inside r, expand. 

An alternative method was given by Squire (2). He assumed that 


the rotational fluid, originally in the region which is ultimately filled by the 
fluid from the source, is completely swept away so that the conditions in 


[Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 2, 1958] 
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the outer region are unaffected by the source flow. Assuming also that the 
total head of the source flow is fixed and that the volume flow from the 
source is the maximum possible, it was then shown that 


r, = (3), Uy = Qr,/V3 = Qr,/v2, (1) 
where Q is the angular velocity of the outer fluid and u, the velocity of the 
inner fluid at infinity. 

These results differ appreciably from Barua’s, which are 


Qr, 


r, = 1-6797r9, Mo = T5098" 





In (2) use is made of the principle that the mass flow attains its maximum 
value, whereas in (1) the mass flow is regarded as fixed. The relation between 
this maximum flow principle and the speed of a long wave propagated 
along the dividing surface is pointed out below, and a result emerges which 
is consistent with equation (1). An account of the calculation of the shape 
of the bulge is given in section 3. 


2. Relation between long waves and maximum mass flow 


The use of the principle that the mass or volume flow takes up its maxi- 
mum value is encountered in many fluid motion problems. Binnie (3) has 











V-w *18% a oy 
> _—_— 
atietlahens AXIS_ OF a 
SYMMETRY 
Fic. 2 


shown that this principle is equivalent to identifying the velocity u, with 
the velocity of propagation of a long wave along the dividing surface, and 
this affords an alternative method of calculation. 

Let V be the velocity of propagation of a long wave of small amplitude 
along the surface dividing a region of solid body rotation from a cylindrical, 
stationary core. Let w be the particle velocity and dr, the displacement due 
to this wave. Then if the waves are brought to rest by superposing a 
velocity V on the inner fluid, as in Fig. 2, continuity shows 

(ro +5ry)*(V—w) = Vr. 
We now equate the total heads at the crest and ‘node’ of the wave to obtain 
bp(V —w)?+ pot $pQ(ro+5ro)? = dpV?+po+4pQ*r5, 


5002 .42 P 
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where p is the density of the fluid and p, refers to the pressure datum in the 
rotating field. Since dr, is small these two equations can be linearized. 
Eliminating w from the resulting equations immediately leads to 

V = Qr,/Vv2 = tp, 
which agrees with (1). 


3. Determination of the boundary between the inner and outer 
streams 


The shape of this boundary was calculated, under Squire’s assumptions, 
by the use of relaxation methods. The problem presented here is that of 
free streamline flow which can be formulated as below. The relaxation pro- 
cedure can be found in (4). Ifq’ is the total velocity just inside the surface of 
discontinuity, then by equating the total head here to that at the rim of the 
bulge, we are led to (q')? = Q2(r?—r,)? (3) 


as the condition to be satisfied on the bulge; r, is the radius of the bulge. 

We introduce the Stokes stream function %’, where since the flow in the 
inner region is irrotational, 

ao , “oe , a , 

o ay = 14 
ca go ie = 0 (4) 
And ! nd ar 

Ox er r or 


where x, r, 6 are cylindrical polar coordinates. Introducing the non- 
dimensional quantities 


r q ub yp’ 


y=-5 (5) 


ry .. Qr,’ 
as he Le 
a4 Lab _ 
0g? " dn? =, On 
throughout the irrotational field, and 


a= HAs (A) = 1 c) 
NBL\ OE ©n} So 


on the dividing surface, where the suffix b refers to conditions there. The 
boundary conditions are 


%’ = %=0 on the axis 


Qr? 
——= OF = — 
3v3 ty 3v3 


=— - 3? 
Qr; 


(4), (3) reduce to (6) 


f, = on the bulge. (8) 


The contour of the bulge obtained is shown in the table below, and Fig. 3 
shows a rough picture of the flow. 


TABLE 1 
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No significant comparison can be made between these results and those 
of Barua’s, since even on the basis of equal volume flow the radius r, and 
the velocity u, would have different values in the two cases. 
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INCOMPRESSIBLE FLOW PAST QUASI-CYLINDRICAL 
BODIES AND SOME ASSOCIATED PROBLEMS 


By L. E. FRAENKEL (Imperial College, London) 
[Received 31 May 1957] 


SUMMARY 
This paper is in two main parts. The first is concerned with axi-symmetric 
incompressible flow past a body of infinite length and of nearly constant radius. 
This problem is reduced to the study of a cylindrical harmonic, G,(x, r), whose 
normal derivative on a cylinder of unit radius is +} (+ for z 2 0). More general 
boundary conditions are satisfied by superposition. In the second part of the paper 
this function is used to derive an approximate solution for the problem of flow past 
slender bodies of revolution the slope of whose meridian section is discontinuous. 


1. Introduction 


THE main object of this paper is to provide a method of calculating the 
incompressible flow past slender bodies of revolution with discontinuous 
profile slope, at zero incidence. This problem seems of interest for two 
reasons. First, bodies like cone-cylinders in transonic, supersonic, and 
hypersonic flow have been studied extensively, and it seems desirable to 
complete the picture by considering the low-speed case. Secondly, the 
incompressible flow past a double-cone provides an approximate model 
for the problem of the vertical water-entry of a cone. The present paper 
is concerned only with the aerodynamic problem; details of the application 
to the water-entry problem are reported elsewhere (1). 

It is clear that in incompressible flow no distribution of sources along the 
axis can provide the required discontinuity in streamline slope away from 
the axis: nor is the body in general a coordinate surface of some system 
leading to separated-variable solutions of Laplace’s equation. Thus classi- 
cal methods of solution cannot be applied. However, the present problem 
can be solved at once by a method developed by Lighthill (2) for the 
corresponding problem in linearized supersonic flow: this method has also 
been extended by Lighthill (3, p. 453), Ward (4, p. 214), and Warren and 
Fraenkel (5). In the present case our first task is to calculate the external 
flow past a quasi-cylindrical body of infinite length, that is, past a body of 
radius R(x) and mean radius ‘a’, such that (1— R(x)/a) is small. The dis- 
turbance velocity potential ¢ can be visualized as due to a source distri- 
bution on the surface of the mean cylinder, and discontinuities in the normal 
derivative now cause no particular difficulty. Our second task is to com- 
bine this solution with that for a smooth, slender body (see, for example, 

{[ Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 2, 1958] 
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4. p. 187) in order to obtain a first approximation to the flow past slender 
bodies of revolution with discontinuous profile slope. 

The solution of the quasi-cylinder problem has other applications. Thus 
¢ can be interpreted as the total velocity potential of a flow emerging from, 
or entering, a porous cylinder with assigned normal velocity; or as the 
temperature in an infinite body with a cylindrical hole, when the steady 
rate of heat transfer at the inner boundary is given. We examine a cylindri- 
cal harmonic which we call G,(z,r), and its partial derivative G,, = G(z,r). 
The function — G, is the response (¢) to a unit step in the normal derivative 
(d,) at « = 0 on the surface of a cylinder of unit radius. The response to 
more general disturbances can then be written in the form of a single super- 
position integral. Asymptotic series are found for G, far from and near to 
the singularity at x = 0, r = 1, and the functions G(z,1), G,(z,1), and 
G(0,r) are tabulated. 


2. The quasi-cylinder problem 
Taking cylindrical coordinates (x,r), consider the problem defined by 


Pent rrt “by = 0, (2.1) 
¢>0 asx>oorr>o, (2.2) 
$,(x,a) = f(x), (2.3) 


where subscripts denote partial differentiation, and it is supposed that 


@ 


| f (x)| dx exists. This represents incompressible flow past a quasi-cylinder 


of infinite length if the total velocity potential is ® = U(x+-¢) and if f(x) 
is the small profile slope. Solutions of (2.1) which vanish for r—> 00 are 
e* K,(\w\r), where A, is the modified Bessel function of the second kind (6). 
Hence the solution of our problem is 


w 


1 ef Ky(\w |r) wir) teed - 
d(2,r) = —2 | Eo tw | f(g) dé (r>a), (2.4) 


> ae 
since Fourier’s integral theorem shows that the boundary condition (2.3) is 
satisfied. Now let 





emtKy(\w|r) 7 3 [ cos wa Ko(wr) dw = G(x,r) (r >1), (2.5) 
— K = ) = ’ = *)s 


se w K,(w) 


—-@®D 


@o 


| oe. dt = eee dw =G,(x,r) (r>1). (2.6) 
1 
0 0 
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Note that 
~6,,(2,1) =} —_ @ = -4 (+ for x 2 0), 
ii 0 
which shows the physical significance of G,. Equation (2.4) may now be 
written ee 


@ 
r («—é r - z—ér - 
$= [o(=*.7\peae-—« [af 7)ye@, es) 
where the last integral is a Stieltjes. 

We proceed to examine the behaviour of G(x,r) and G,(x,r) in various 
parts of the field. For aerodynamic applications the values on the body 
(r = 1) are probably the most important, but in the water-entry problem 
G(0,r) (which contributes to the normal velocity on the free surface) is 
also of interest. 


3. The cylindrical harmonics G and G, 

3.1. Asymptotic expansion of G, for x*+-r® + oo. 

Before calculating the asymptotic behaviour of G, far from the origin, 
we evaluate three relevant integrals. These are all based on the well-known 
result (6, p. 388): 


if ; ‘Se, ee 
me ! COS WX K,(wr) dw = 2 (22-7) == 2p’ (3.1) 


0 
(where p? = x*+-r? in this section), and are 


A,= sin wr wl Ky(wr) dw = } sinh! 


l 
7 


, x 
sin wx w Ky(wr) dw = — 





> , x 2p? e 
sin wx w log w K,(wr) dw = 3 log-—+y—2 


where y is Euler’s constant. A, and A, follow from integration and differen- 
tiation of (3.1) with respect to x; the third integral is treated in Appendix I. 

Expansion of 1/K,(w) in the integrand of (2.6), in ascending powers of w 
and log w, now gives 


G,(x,r) = A,-+4log 2—y+)4q—}4g+-O(p~log*p) 


r, 
- 


~~ « 
= jsinh-1= 4 7 “(log “E—5\+0(0 log?p). (3.3) 
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In what follows we examine @ rather than G,, because this happens to 
be more convenient in subsequent applications: G, can be recovered by 
simple integrations. 


3.2. Asymptotic expansion of G forx>0,r->1 

To calculate the behaviour of G(z,r) near the singularity at (0,1), it is 
convenient to introduce the analogue of G in the supersonic problem. We 
write x—r+1=t, r—1 = y, so that ¢ is constant on the characteristic 
surfaces of the supersonic problem, and introduce a function V(t, y) defined 
by 





> 1 “T ept+K.(p-+-py) 
Vit.y) = 5h | Rp? (>) (3.4) 


eK (p+Py) _ » | eV (t,y) dt (rep > 0). (3.5) 
K,(p) : 
The properties of V(¢,0) = U(t) are well known: the function and its deriva- 
tive are tabulated (7); for small ¢ (4, p. 173) 

U(t) = 1—}t+ C—O + Bhttt+...; (3.6) 
and for large t, U(t) ~ t-!. Expansion of the integrand in (3.4) in ascending 
powers of y gives 

Vit,y) = Ult)+yU'(H)+y¥U"(H+y[FU"(O—jU"()] +... (3.7) 
dashes or superscripts denoting differentiation. Returning now to G, we 


use (3.5): a 
ie 1 [ e*K,(wr) 


G = re—  jeaieameaianeeae 
(x, r) rex | aK, (a) dw 


= re ett day e-ev | eV (t, y) dt 
0 


0 
a= #02 | VGY) ay, (3.8) 
7 t—iz 
0 


where z = x+1y. 
In view of (3.7) the problem is now that of expanding the integral 


U™(t) 


eM (3.9) 


B(z;m) = | 
0 


in ascending powers of z and logz. Writing 

fof 2 1 1 

| .. | logz (dz)” = = |logs—(14+5+--+2)| = L,(z), 
0 


n! 
0 
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we obtain by repeated integration by parts 


N ° 
a re 1y"41Z, (—iz)Um+0) +-(— 1441 Ly(t—iz) UN +t) dt. 
_— 0 


Here the series is of the required form. The last (integral) term, E(z;m, NV) 
say, has N derivatives with respect to z at z = 0: 


B(0;m, N) = (—1)8(—i)" [ Ly_,(t}U™+*W) at 


0 


=. 


= —i" [loge Uomm+n(t) dt 
0 
= 1"... Say. (3.10) 


Hence expanding £ in a finite Taylor series, and then letting N — 00, we 
obtain 


B(z;m) = > | (—1)"#1LZ,(—iz)U™+""(0) + Cn oe (3.11) 
at n! 


Note that the logarithmically singular part of B is contained in the first 
group of this result, the values of U’™+™(0) being known from (3.6) or 
(3.4). The second group is regular, and the coefficients C,,,,, have to be 
found numerically. The first two required here are 
C, = 0-4568, OC, = —0-164. 
It now follows from (3.7) and (3.8) that 


G(a,1) = I og a(—1 + fx? — Gh0*)+C,+-4rx— 
7 


—($+4C,)a?@—Jdnx3+O(2x4)| (x > 0) 
G,(x,1) = * [log 2(—a-+pya®—,fhyr5)+(1+-0,)2-+ 
+ §ra*— (§+3C,)a*—hera*+ O(a*)] (x > 0), 
G(0, 1+y) [log y(—1—fey?+ dey*) +Cy+ 4y—(&—4C)y? + 
+(&—IC,)y3+ Ol(ytlogy)] (y > 0). 
3.3. Numerical integration of G(x,1) and G(0,r) 


«x 


The integral G(x,r) = re 1 [_ Pm 
mJ wK,(w) 
0 


while convenient for x = 0, r > 1, is oscillatory elsewhere and does not 
converge rapidly on r = 1. Accordingly, we take a contour, in a complex 
w-plane, consisting of the positive real axis, the positive imaginary axis, 
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a small quarter-circle at the origin, and a large quarter-circle at infinity. 
This contour contains no singularities, and the contributions of both 
quarter-circles vanish in the limit, for z > 0. Thus one finds that 





(3.15 a) 


vy 1 fe Spr) ¥Q(A)—J(A)KQ(Ar) 
G(xz,r) = = (< = FRA) LYRA) da, 
0 


(3.15 b) 


Ge.) == [5 TFT 


Note that G(x,r), which is defined only for r > 1 by (2.5), is continued 
analytically by (3.15a) into the domain xz > 0, r > 0, the point x = 0, 
r = 1 being excluded; G, may be continued similarly. 

Equation (3.15 b) has been used for the numerical evaluation of G(x, 1). 
If asymptotic series for the Bessel functions are substituted for A > L say, 
the integral from L to oo can be written as a rapidly convergent series. 
G,(z, 1) has been found by numerical integration of G. These results are 
shown in Table 1, and G(0,r) is shown in Table 2. In addition we have the 
series (3.12) to (3.14) for small values of the argument, and the following 
series (cf. equation (3.3)) for large values: 


G(x,1) = ——3 (log 2x—1)+O(x-log*x) (x > 0), (3.16) 


G,(x, 1) = 4 log + 2 (Hog 2x—4)+O(a-‘log*x) (x > 0), (3.17) 


G(0,r) = 1, gain (3.18) 


G is an even, G, an odd function of x. A comparison of values found 
by numerical integration with those given by the various series is shown 
in Table 3. 


4. Slender bodies of revolution with discontinuous profile slope 

Let R(x) be the radius of a slender body of revolution at zero incidence, 
of length / and in a uniform stream of velocity U; let the cross-sectional 
area be S(x) = 7R*(x). The ends of the body are either pointed or cylindri- 
cal and extending to infinity. The length / is assumed to be of O(1), and 
R, R’, R" are of O(e) everywhere, except that at the points x = c, there 
are discontinuities in R’ and/or R” such that 


S’(c,+)—S’(e,—) = AS;, S’(c,+)—S"(e,—) = AS}. 
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TABLE TABLE 2 








G(x, 1) G(x, 1) G(o, r) 


™ 








co ° 
0°9013 01208 
o-7Oll 0°1997 
0°5907 0°2639 
05163 o"3190 
04616 0°3678 
04188 O'4117 . O13Is 
Oo: o-4518 . 0°0853 
0°3555 04888 . 0°0634 
0°3: o°5231 0°0503 
o3 o°555! 


@ 

0°4477 
0*3267 
0°2447 
o18i1 


ee eee 
oumouwmo 








O°5551 
06135 
06659 
07133 
07566 
07966 
0°8337 
08684 
09008 TABLE 3 
o°9313 
00-9602 


2°8 
3° 


e2929290929°9909 





Function | Numerical integration Series 





3°0 , 09602 G(o-4, 1) 05163 05163 
4°0 . 10841 G(o-6, 1) 0°4188 o°4179 
50 , 11840 G(10°0, 1) 0°0482 20480 
60 "07 1°2677 
70 , 1°3397 G,(10°0, 1) “5095 15103 
8-0 "0! 1°4028 
g°0 "053 1°4590 G(o, 1°5) 7 0°4496 
10° 048 1°5095 G(o, 10°0) 00505 




















Note: Entries for z = 0 to z = 0-4 
are based on the series (3.12), (3.13). 


S'(x) in this problem corresponds to 27af(x) in the quasi-cylinder 
problem: hence by analogy with (2.8) we postulate the disturbance potential 


he 1 ry a—é fr - 
d(x,7) = — | aS me) dS'(é). (4.1) 


This certainly satisfies Laplace’s equation, and in Appendix II it is shown 
that ¢ satisfies the boundary condition 


$,(x, R(x)) ~ Riz), 


the error decreasing continuously from that of linearized wing theory near 
a discontinuity, to that of smooth-slender-body theory away from it. 
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By (4.1) the axial disturbance velocity is 


Mc r za—é fr dS'(£) 4 
Pale,9) = =} OCR Re) RO R€)’ eat 


—@ 


but it follows from arguments like those of Appendix II that this may be 


Pal) ~ sim J ore Fal yt a 


without further loss of accuracy. Then on the body 
— As.+ ¥ @,[%—%,1)as7 
¢.(x, R(z)) sl Bah 2 “Ra 1) i > Rey )a ms, 


4 { 0, 1)s"@) ae], (4.4) 


replaced by 


and in the last integral of (4.4) one may substitute the asymptotic form 
(equation (3.17)) G,(x, 1) ~ +4log|2z|. These results are of course wholly 
analogous to those of the supersonic theory (2, 3, 4, 5). 

If there is a discontinuity in slope (and in S’) at x = 6, R(0) > 0, say 
AR'(0) = e, then by (4.4) and (3.12) the total velocity on the body is locally 


i ~ 1+ <log|x|+O(e loge). 
j 7 


If the pom nose of the body is at = 0, R(0) = 0, say AR’(0) = «, 
AS"(0) = 2me®, then by (4.4) and (3.17) the total velocity on the body is 
locallyt 
i ~ 1+} log z+ O(c? log «). 

As Van Dyke (8) has pointed out, such results are formal expansions (for 
small «) of the singularities which occur in exact theory. The exact solutions 


are locally 
i~ |ar |e +O%e*) and i~ atet+ Ofetlog «), 


The present results can be made uniformly valid by Van Dyke’s method, 
but as the regions of invalidity are exponentially small this is probably 
not worth while in practice. 


+ This value does not come from the AS*(0) term in (4.4), but from the rest of 


1 


“hs Cre 1) ase, 


in which R(x) ~ ex + 0, G, ~ —}log[(é—z)/ez). 
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As an example of the method of this section we calculate the pressure 
distribution on a cone-cylinder, defined by 
Rixw)=ex (0<2x<1) 


=e (i < @). 
(4.4) becomes 


] xr—1 x xr—l1 
= ¢*/—__. G[ —__. , 1} — G,|—_. , 1] +. G, |—_., ] } I, 
nee = Fours ) OR) }+ Fee ) 
and the pressure coefficient is 
a] P—Po 5 
Co = t29 = — — 2, 
Pp bpU? 2¢, ¢; 


where ¢, = R(x) on the body. Fig. 1 shows the pressure distribution for 


e= 01. 
1-0 6 s 24 
€ 
z 


Late 
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APPENDIX I 
THE INTEGRAL A, 


To evaluate the integral 


A; = | sin wr w logw Ko(wr) dw, 
0 





INCOMPRESSIBLE FLOW 
consider ‘ | sin wr w~* Ky(wr) dw = M(x,r,k), say (k < 2). 
0 
Clearly Ay = Myz;(2,7, 1) 


(subscripts denoting partial differentiation). By Fourier inversion of the basic 
relation (3.1) 


M(x,r,k) Li 
7 


sinasortde { cos wt (€2 + r*)-t dé 


0 


> (weg eel 
J) +r sin $k (k—1)! 


(+ fora 2£,0<k<2,2> 0). 


va log(x+£)-+log|x—€| (ae 
Hence M,{ z,T, 1) = aa if (€2+r)t dg tty | (+r) 


Care must be taken in differentiating this expression with respect to 2, because of 
the singularity at € = x: when this is done M,, may be found by elementary integra- 
tion, and a further differentiation yields A, (equation (3.2)). 








APPENDIX II 


THE BOUNDARY CONDITION (4.2) 


In this appendix we examine the boundary condition satisfied by 


_1 fg (ee 
ser = —3 | Olea RG) oe yeti 


-@ 


Consider first the contribution of a single discontinuity, which we take at x = 0, 
R(0) = R, for convenience: 


bar) = —F alz. z). 


Two regions are studied separately: 0 < |x| < O(e) and O(e) < |x| < O(1). 
(i) 0 < |x| < Ole). 


AS; 
By (2.7) dolx, Ry) = +o 
0 


and since R(x) = R,[1+ O(x)], while do,, is O(ex-!) at worst, 


(+ for x 2 0), 


MN AS; 
dor(z, R(x)) = + GrRa) tt +O). 
(ii) Oe) < |x| < O()). 
With x large, r of O(1), (3.3) becomes 
2\z| 


G,(z,r) = +[dlog7#! 4 2 (5 +10g 4-3) +...], 
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so that here 


AS; 2 R3 (1 a 3 
o(2.r) = F = [410g 145 5 * Ratios; R 3) + 


427\2 


1 1 RR} (R(x) 
2 Ra) (te) sing 


_ AS5 
Por( =, R(x)) = F —|- 
and since R(x) = R,[1+O0(z2)], 
ASo 
x, R(x 
(If the body length is greater than O(1), e.g. if the body becomes cylindrical asymp- 
totically, then the error factor for x > O(1) is 1+ O(e*z~?).) 
Now ¢ is a sum of finite and infinitesimal terms like ¢,: to generalize the analysis 
above we need merely replace AS} by dS’(€) and sum. The intervals (i) anu (ii) above 
now refer to x--& instead of x. Then 


pki + Ole 1) ]). (II. 3) 


5 ] p- de j . 
$,(2, R(x)) ima | - dS’(€)[1 + O(e%(x—£)) 


es 
r 2+ Ae 


+ f ~ ase , (II. 4) 


2~—Ae z 


2S'(xr) 
i 4n R(x) 


R’(x) 


where A is some finite number. The error associated with the discontinuity terms 
is given by (II. 2) and (II. 3), with 2 replaced by x—c;, and by (II. 4) the error asso- 
ciated with S*(&) is a factor 1+ O(e* loge). 

Finally, we remark that in going from (4.3a) to (4.3b) it is legitimate to replace 
R(é) by R(x) because for (w—£)/R(E) < O(1), R(é) ~ R(x); while for 


(x—£€)/R(E) > O(1), 


the leading term of wie (a ae) 


is independent of the reference length R(€). 
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ON THE CALCULATION OF EDDY VISCOSITY AND OF 
HEAT TRANSFER IN A TURBULENT BOUNDARY 
LAYER ON A FLAT SURFACE 


By D. E. BOURNE and D. R. DAVIES (University of Sheffield) 
{Received 25 April 1957] 


SUMMARY 

This paper describes a method of calculating the distribution of eddy viscosity 
and of rate of heat transfer in forced convection from a flat plate through a develop- 
ing, incompressible, turbulent boundary layer. It is necessary to know (a) the main 
stream velocity and plate temperature (both assumed to be independent of down- 
stream position), and (6) empirical formulae for the distribution of surface shearing 
stress, temperature in the laminar and transition sub-layers (involving the Prandtl 
number a), and mean velocity in the fully turbulent, transition, and laminar layers. 
A similarity law for the distribution of mean velocity in the fully turbulent layer 
is calculated from these empirical formulae, and is used in an integration of the 
boundary layer equations of mean flow. This enables the distribution of Reynolds 
shearing stress and of eddy viscosity to be calculated. The distribution of eddy heat 
diffusivity is then available, using Reynolds analogy, and the method of calculation 
developed recently by Davies and Bourne (1) is applied to evaluate rates of heat 
transfer in the case of air (o = 0-72). These calculated values are found to be in 
reasonable agreement with measured values given by Elias (4). Finally, it is suggested 
that this new method of calculation may be applicable in other important problems, 
involving convective heat transfer, where measured velocity profiles are available 
and eddy viscosity can be computed from the equations of mean flow. 


1. Introduction 


ReEcENTLY Davies and Bourne (1' have shown how the distribution of rate 
of heat transfer from a uniformly heated flat plate into a developing turbulent 
boundary layer can be calculated, when the plate temperature is indepen- 
dent of downstream distance. The measured mean velocity profile over the 
plate is assumed to be available; the distribution of eddy viscosity is also 
required, and can be deduced from values of eddy viscosity measured by 
Townsend (2) over a smooth wind tunnel wall. 

It is important to note that in this calculation (1) three assumptions are 
made. 

(i) Itis assumed that Reynolds’ analogy .s true. According to this assump- 
tion, the transfer of momentum in turbulent flow is exactly analogous to the 
transfer of heat and matter. In the case of heat transfer, ifr, is the Reynolds 
shearing stress at a normal distance y from the plate, q, is the rate of heat 
transfer at the same point in a direction perpendicular to the plate, p is the 
fluid density, c,, the specific heat, 7’ the mean temperature, and U is the 

(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 2, 1958] 
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mean velocity parallel to the plate, then Reynolds’ analogy may be expressed 
by the relation 
p(eUjey) pe, (eT /ey) 

where « is the eddy viscosity. Reynolds’ analogy has been experimentally 
verified for forced convection in turbulent flow near a solid surface (see, for 


example, a discussion by Davies (3)). 





(ii) A particular similarity law is assumed to exist for the distribution of 
mean velocity in the fully turbulent part of the boundary layer, i.e. the 
velocity is assumed to be a function of a similarity variable é of the form 


3600 











-_ 
Uv 
7) 
“ 

~~ 
“ 
E 
7] 
— 
~_ 
~ 




















00 
— O01 £ 002 


Fic. 1. Mean flow distribution U (em/sec), measured by Elias 
(U, = 35 m/sec) and calculated from empirical logarithmic laws. 
—— mean curve through the measured values (U = 1-72U, £%); 
---- calculated from Coles’ logarithmic formula, at x = 20 cm; 
cocnence calculated from Coles’ logarithmic formula, at z = 50 cm 
y/a4, where x is the downstream distance from the virtual origin of the 
turbulent layer and q is an experimentally determined empirical constant. 
Experimental results presented by Townsend (2) and an analysis by Davies 
(3, Fig. 1) of experimental results by Elias (4) show that this assumption 
is valid over the limited range of x values considered. 
(iii) The validity is assumed of a set of empirical formulae for the distri- 
bution of mean velocity in the laminar, transition, and fully turbulent zones 
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of flow near a flat surface. These formulae are presented and discussed by 
Howarth (5, p. 824) and more recently by Coles (6). The derivation of the 
corresponding formulae for the distribution of temperature in the laminar 
and transition zones is discussed also by Howarth (5, p. 826); the shearing 
stress and vertical heat flux is assumed to be independent of y in the sub- 
layers and the temperature is expressea as a function of Prandtl number co, 
kinematic viscosity v, surface shearing stress 7), and height y above the 
surface. 

This method of calculation is now extended in order to make predictions of 
heat transfer rates, when measured values of the mean velocity profile and 
of eddy viscosity are not available. This can be done by first applying a 
well-known empirical formula, given by Goldstein (7, p. 362), to calculate 
the surface shearing stress in terms of U,, the main stream velocity, and then 
using the logarithmic mean velocity formula for the so-called ‘law of the 
wall’ as given, for example, by Coles (6). Using these formulae, the velocity 
distribution over the range of downstream distances to be considered is 
calculated and plotted as a function of the similarity variable £. In the 
conditions investigated by Elias (4), involving forced convection from a 
flat plate 50 cm in length, good approximate similarity of the calculated 
velocities is displayed in the inner 30 per cent of the boundary layer, if a 
suitable value of q is chosen. The calculation of mean velocity distribution 
in the outer parts of the layer could be improved by using the so-called 
‘wake function’ suggested by Coles (6); this calculation would be compli- 
cated by the fact that the wake function depends on the degree of turbulence 
in the oncoming stream and on the ‘intermittency’ in the outer part of the 
boundary layer. However, it is the distribution of flow variables in the 
inner 30 per cent or so of the layer which is crucial for heat transfer calcula- 
tions and, as shown by Davies and Bourne (1), the distribution of flow 
variables in the outer parts does not appear to affect significantly heat 
transfer at the surface of the plate. The application of the law of the wake 
has consequently not been considered in this paper; indeed, in the con- 
ditions of Elias’s measurements this is found to be a comparatively small 
effect and in this case can be neglected for calculations involving only the 
inner part of the layer. 

The calculated similarity velocity profile in the inner part of the boundary 
layer over the flat plate is found to correspond closely to a mean profile 
drawn through the velocity measurements made by Elias and is used in the 
calculation of eddy viscosity «. Following a similar idea introduced pre- 
viously by G. I. Taylor (8) in calculating eddy viscosity and mass transfer 
in turbulent flow through a pipe, the Reynolds stress distribution and hence 


the eddy viscosity distribution over the whole field of flow is calculated by 
5002 .42 Q 
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integrating the equation of mean flow using (i) a numerical method based 
on the actual similarity velocity profile, and (ii) an alternative method based 
on a power law representation of this similarity profile. The eddy viscosity 
distribution obtained is found to be in reasonable agreement in the inner 
part of the layer with that previously deduced by Davies and Bourne (1) 
from measured values given by Townsend (2), and the ensuing calculated 
values of rate of heat transfer are also in reasonable agreement with meas- 
ured values given by Elias (4). 


2. The approximate similarity profile in the inner part of a turbu- 
lent boundary layer 
Measurements of mean flow distributions in turbulent boundary layers 
have been discussed by many investigators (see, for example, Coles (6)) 
who show that the logarithmic formula for velocity (parallel to the surface) 
U/U, = A log,(yU,/v)+ B (1) 
gives a good representation in the region of applicability near to the wall, 
which is the region of greatest importance for calculating heat transfer 
from the wall; U, = ,/(79/p) is the friction velocity, v is the kinematic 
viscosity, and A and B are constants. Slightly different values of A and B 
have been suggested by Coles and Clauser for flat plate flow in the absence 
of mean pressure gradient and by Nikuradse for pipe flow, and are discussed 
by Preston (9). We find, however, that calculated values of heat transfer 


based on these alternatives are all within reasonable range of measured 
values by Elias (4). 

An empirical formula giving U, in terms of U}, v, and x has been developed 
from results of pipe flow [see, for example, Goldstein (7), p. 362], and is 
expressed in the form ' 


U2/U} = 0-0289[v/(U, x)}*. (2) 
At the highest main stream velocity considered by Elias in his flat plate 
measurements (U, = 35 m/sec; 1, length of plate, = 50 cm), an analysis by 
Davies (3) of the velocity values suggested that the completely laminar 
layer near the leading edge of the plate was not appreciable and the virtual 
origin of the turbulent layer could be taken along the leading edge. A 
numerical evaluation of U,(x) in these conditions was carried out by two 
methods, firstly using Coles’ form of equation (1), i.e. with A = 5-75, 
B = 5-1, in conjunction with the actual measurements of U by Elias, and 
secondly using empirical formula (2). The values of U, given by the first 
method are about 3 per cent larger than those given by the second method, 
although the dependence of U_ on x is about the same; the first and second 
methods give U, to be proportional to z-°™ and z~°” respectively. The 
same magnitudes of U’, are given, to a sufficiently good approximation, by 
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both methods over the range of x considered, if the proportionality factor 
in pipe flow, 0-0289, is altered to 0-0300—a somewhat similar correction to 
that suggested by Goldstein (7, p. 362). This suggests that the empirical 


fi Y y € / T . 
orm U2/U? = 0-0300[v/(U, x)}* (3) 


is a good working formula to use over flat surfaces in order to compute U,, 
provided U’, is large enough (and / long enough) to neglect the completely 
laminar zone near the leading edge; we note also that the validity of (3) 
appears to be further restricted, according to Goldstein, to Reynolds 
numbers of order 10*, a condition which is satisfied in the case to be 
considered. 

Taking x = 20, 30, 40, 50 em with U, = 35 m/sec and with appropriate 
values of U,, the velocity values U given by the logarithmic formula (1), 
with A = 5-75and B = 5-1, were then computed. These values are plotted 


against a suitable similarity variable, which can be determined from 
equation (3) and the analysis of the heated plate problem, given previously 
by Davies and Bourne (1, equations 24-31). On the basis of this analysis 
we see that the friction temperature 7’, is independent of z and, since equa- 
tion (3) shows that U, is proportional to 2~-®?*, it follows that the local flux 
of heat, given by q, = pe, U,T, (where p and c, denote the appropriate 
density and specific heat), is proportional to z-°!®, The analysis (1, equation 


31) also shows that q, is proportional to 2*-!, where q is the index in the 
similarity variable '€ = y/(ket). (4) 


It follows that we must take g = 0-90, which is in close agreement with the 
value 0-89 derived previously by Davies (3) from an analysis of the velocity 
measurements of Elias (4), and leads to a variation of g, with z in close 
agreement with the measured variation reported by Elias. Taking q = 0-90 
and k = 1-0, the values of U, when plotted against £ for various x values : 
(shown in Fig. 1), demonstrate reasonably good approximate similarity for 
values of x between 20 cm and 50 ecm in the inner 30 per cent of the layer. 
A mean curve drawn between the calculated points corresponding to 
x = 20cm and 2 = 50 em fits very closely the power law U = 1-72U, €°"; 
we note that the same mean profile is obtained if we use the values of A 
and B suggested by Clauser. This power law also fits quite closely the 
velocity profiles U(€) measured by Elias and by Townsend (q = 1, k = 1), 
and hence we can follow our previous method (1) and deduce the distribu- 
tion of eddy viscosity appropriate to a flat plate from the measured results 
given by Townsend for a wind tunnel wall: this distribution can then be 


compared with the numerical values calculated from the equation of 
mean flow. 
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For values of x greater than 100 cm, the values of U computed from 
equation (1) fall below the similarity curve shown in Fig. 1 for all values of 
€, and this clearly indicates that the method of calculation suggested in 
this paper must be restricted to a range of Reynolds numbers which can be 
determined only by computing U from equation (1) and plotting against 
the similarity variable é in the inner 30 per cent, say, of the layer thickness, 

As seen in Fig. 1, the wake function introduced by Coles (6), which is 
essentially the difference between the measured values of U and those 
given by the logarithmic law of the wall, is small even in the outer part of 
the boundary layer and plays no part at all in the inner zone. We note, 
of course, that at larger downstream distances it may become important 
and may have to be included in heat transfer calculations. 

The similarity profile, shown in Fig. 1, is now used to calculate the 
distribution of eddy viscosity from the equation of mean flow. 


3. Calculation of Reynolds stress and eddy viscosity distribution 

Let (U,V) and (u,v) denote components of mean velocity and eddy 
velocity respectively in a fully turbulent boundary layer, relative to co- 
ordinate axes Oz, in the plane of the plate and in the direction of the main 
stream, and Oy, normal to the plate. The equation of mean flow, when the 
pressure gradient is zero, is 


,oU_ ,,0U RP nee - 
U—+4+V— = —-< (pa), (5) 
dx oy poy 
where p is the fluid density and pz the Reynolds shearing stress; we neglect 
p(eU /dy) in comparison with pv? in the fully turbulent zone. In the incom- 


pressible case the equation of continuity of the mean flow is 


UW» 
Ox cy 
so that we may define 
Gin wt Wa 2. (6) 
cy ba 


where (x,y) is the usual stream function associated with the mean flow. 
We now write 

ob = xf (€), (7) 
where £ is the similarity variable defined by equation (4). Equations (4), (5), 
(6), and (7) then lead to a simple differential equation for (7), 


o 


ag) pan gk-txa-p EF (8) 





—~ #f- mm © 2 
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The function f(€) is evaluated from the mean similarity profile computed 
in section 2 of this paper. We have 


U=* = ir-y¥'o) (9) 


primes denoting differentiation with respect to £, and the function f’(€) is 
thus known numerically. The derivative 


f’ = kU’ (10) 
may be computed by numerical differentiation and the integral 


é 
fak [Ud (11) 
é. 


by numerical integration; €, refers to the assumed lower limit of the fully 
turbulent zone, following the method of analysis applied by Davies and 
Bourne (1). Substituting (10) and (11) into (8), and integrating (8), we find 


é é 
i = qkxt- [ uf U ae) dé +C, (12) 
é : 
where, following Davies and Bourne again, C is evaluated from conditions 
on the upper boundary of the transition layer at the downstream distance 
x” == x, say, at which heat transfer is being considered. The height, y,, of 
the transition zone at x = x, can be calculated from the empirical logarith- 
mic velocity form. If this point corresponds to = €, = y,/(kx#), say, and 
if the associated value of (@) is (@7),, then 


€ é 
(i) = (WH), + ghee? f U'{ f UV ae) ae. (13) 
fs fs 
. eu ef? 
From equation (9) — = k-*z-af", (14) 
cy 


and consequently the eddy viscosity «, defined by —(@)/(@U /ey), may be 
expressed in the form 


£ 


£ 
€ _ K@),2'* gh ee ; ms 
U, ah2a-) ra ee x ea = uF U U dé) dé. (15) 


At a point just above the upper boundary of the transition layer the 
contribution to the total stress due to the viscous stress can be neglected, 
and it can be shown that the variation of total stress in the sub-layers can 
also be neglected, so that to a sufficiently good approximation 


(Wd) » = —- = — U3, 


where 7, is the total shearing stress at the surface of the plate. In order 











230 D. E. BOURNE AND D. R. DAVIES 


to verify this constancy of 7, in the sub-layers, we have computed the 
expression for variation of 7, in the conditions of Elias’s experiments, when 
the pressure gradient is zero, given by Coles (10) 


yU-ly 


ee U "a ee] 16 
To AU, | (7) vf mn 
0 
7 T T 
whereA = — 7 *. Using the laminar sub-layer velocity form A == yU, 
T x rT Vv 
and equation (1) with A = 11-0 and B = —2-69 for the transition layer, 


integration of (16) showed that the variation of + 
transition layer can certainly be neglected. 

It has not been found possible to treat the heat transfer problem using 
the form (15) for «, since the right-hand side depends not only on é but 
also on x. However, the variation with x is found to be very slow and can 
be neglected to a good approximation: the calculated values of ¢/U, 4-1 
along the centre (€ = 0-005) of the inner crucial zone are shown in Table 1 
and illustrate this slow variation with 2 over 60 per cent of the length of 
the plate lying upstream from x, = 50cm. In calculating the rate of heat 
transfer, dp, at x, = 50 cm in Elias’s experiments, the values of (#7), 2'~¢ 
and é, at this distance are computed and assumed to be constant on the 
right-hand side of equation (15). The error involved, due to this factor, 
in computing q, at this distance is then likely to be small, as the calculated 
value at x = x, probably depends critically on the distribution of « near 
x, but is not so dependent on conditions further upstream, the main contri- 
bution to the temperature near the surface coming from the heat released 
in the immediate upstream neighbourhood of z,. We find that the ensuing 
calculated values are in reasonable agreement with measured values. 


up to the top of the 


v 


TABLE | 
Calculated values of eddy viscosity along the centre of the inner crucial zone 
(€ = 0-005) in the conditions of Elias’s experiments 


(U, = 35 m/sec, 2, = 50 cm, q = 0-90) 





x (cm) 20 30 40 5° 

















10%e/U, x*a-1 1°24 1°20 V1] erg 





This numerical method of calculation (using equation 15), involving two 
numerical integrations and one numerical differentiation, is considerably 
simplified in cases where a power law of the form 

U/U, = bé8, (17) 


where 6 and £ are constants, can be constructed to fit with reasonable 








ON THE CALCULATION OF EDDY VISCOSITY 231 


accuracy the similarity velocity profile obtained by the method of section 1. 
This can be carried out in the conditions of Elias’s measurements and the 
following analysis applied. Using (17), equations (6) and (7) give 
f = kb(1+8)7U, 6, (18) 

taking ¢ = 0 at € = 0, and integration of (8) with this form for f yields 

(7) = (TB) 4+ kb*gB(1+-B)-"(1 +28) Uf aa EP— EEF), (19) 
so that (%@) = (W@), at € = &,. The eddy viscosity is then expressed in 
the form 


ee. -|-- k*bq 
bpU} (1-+8)(1+ 28) 





(estas) e-#, (20) 
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Fic. 2. Calculated distribution of ¢/(U,x”) in the conditions 
of Elias’s experiments; U, = 35 m/sec, x, = 50em,q = 0-90. 
calculated values, based on the power law method, using 
Clauser’s logarithmic velocity law; ---- calculated values, 
based on the numerical method, using Clauser’s logarithmic 
velocity law; ----- calculated values, based on the power law 
method, using Nikuradse’s logarithmic velocity law; « values 
inferred from Townsend’s measurements 
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The expressions (15) and (20) for «/(U, 2”) have been evaluated, using 
the conditions associated with Elias’s experiments (U, = 35 m/sec, 
x, = 50cm) and Clauser’s empirical logarithmic law for the velocity profile 
(A = 5-6, B = 4-9 in equation 1), and the results are shown in Fig. 2, 
together with values deduced from Townsend’s measurements by the 
method discussed by Davies (3). The numerical values given by (20), when 
based on Nikuradse’s profile for pipe flow, i.e. with A = 5-8, B = 5-5 in 
(1), are also shown; the values (not shown) based on Coles’s recommended 
profile, with A = 5-75, B = 5-1, lie approximately midway between the 
values based on Clauser’s and Nikuradse’s empirical profiles. These results 
are all in reasonably close agreement with each other in the crucial inner 
zone of the boundary layer. The computed values derived from (15) diverge 
markedly from the other values in the outer 50 per cent of the layer, due 
to difficulty in evaluating numerically the slope of the mean velocity profile 
in this region. 


4. Calculation of rates of heat transfer 

In order to compute by a reasonably convenient method (see Davies and 
Bourne, 1) the rate of heat transfer by forced convection (neglecting aero- 
dynamic heating) from a heated smooth flat plate, when the main stream 
velocity is sufficiently large to neglect the effect of the leading edge laminar 
layer, power laws of the form 


€/(U, x?) = ag (21) 


were fitted to the various computed « distributions in the inner 30 per cent 
of the thickness of the layer. As in a previous case discussed by Davies 
and Bourne (1, Fig. 1), this was carried out with considerable accuracy in 
this region. Reynolds’ analogy can then be applied so that eddy heat 
diffusivity ¢«,, = «, and the formulae given by Davies and Bourne (1, 
equations 31-35) used to evaluate rate of heat transfer Q for various lengths 
of plate up to the maximum value of 50 cm considered by Elias. The 
numerical results for Q, based on all the « distributions shown in Fig. 2, 
were found to be within 8 per cent of the measured values given by Elias, 
the extreme calculated values for z, = 50 cm being 0-59 and 0-52 watts 
per degree C (corresponding to the Clauser and Nikuradse profiles respec- 
tively) compared with the measured value of 0-56 watts per degree C. 
We note also that the temperature distribution and rate of heat transfer 
are influenced by the important factors (a) the Prandtl number o, and 
(b) the presence of the laminar sub-layer and transition layer, and their 
effect can be evaluated by the equations obtained in (3) (equations 14, 27, 
28), or by the alternative forms of (1). For example, when o = 1, we find 
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by numerical evaluation of the equations, that the temperature and velocity 
curves are identical, but when o = 0-72 (as in air) they differ appreciably 
(see (3), fig. 2, p. 335). 


5. Conclusions 


The method of calculating the distribution of rate of heat transfer on a 
heated flat plate presented previously by Davies and Bourne (1) has been 
extended in this paper, so that the initial data required are (a) the main 
stream velocity, (b) the surface temperature of the heated plate, and 
(c) empirical formulae for surface shearing stress, temperature in the sub- 
layers, and mean velocity in the fully turbulent layer and sub-layers. This 
new method eliminates the need to know the measured velocity profile and 
the measured distribution of eddy viscosity; the former can be calculated 
from the empirical logarithmic formulae and the latter then evaluated from 
the computed similarity velocity profile. The computed values of rate of 
heat transfer are in reasonable agreement with available measured values 
on a heated flat plate. We note, however, that the method is only applicable 
when U, is large enough to reduce the extent of the completely laminar 
layer near the leading edge sufficiently to neglect it, so that the virtual 
origin of the fully turbulent layer is not downstream of the position of the 
leading edge of the plate. We note also that the application of the method 
must be restricted to a range of x-values for which the downstream Reynolds 
number is of order 10°, so that the empirical formulae may be applied, and 
that the velocity profile must show reasonable approximate similarity, at 
least near to the plate, in this range. 

Finally, it is suggested that this new method of calculation may be 
applicable in other problems involving convective heat transfer (e.g. in the 
turbulent region of the boundary layer over a heated rotating disk), where 
the distribution of eddy viscosities is computable from the available 
measured velocity profiles and the equations of mean flow. 
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SUMMARY 


In this paper a system of transmission lines coupled by reactive networks at regular 
intervals is examined. If the coupling networks are identical, it is shown that modes 
of propagation, equal in number to the number of lines, can exist. These modes are 
typified by voltage and current distributions which can be determined from the 
properties of the network, and show similarities to waveguide modes. In particular, 
the modes are orthogonal in that the power flow when several modes exist is given 
by the contributions of the modes considered one at a time. 


1. Introduction 


PROBLEMS involving two or more transmission lines coupled together at 
regular intervals arise frequently in microwave engineering. An obvious 
example is the directional coupler (1) which in one form consists of two 
waveguides coupled together by a succession of equi-spaced slots. Con- 
siderable interest is now centred on the use of the H), mode in circular 
guide (2), since this has the property of decreasing attenuation with 
increasing frequency. Guides sufficiently large to allow this mode to 
propagate also allow a number of other modes to propagate, and any 
discontinuity will result in coupling between the modes. One method of 
reducing the undesirable effects associated with these other modes is to 
insert a number of equally spaced identical obstacles chosen in a way 
which leads to the suppression of at least one of the unwanted modes. In 
the above two examples, the transmission lines corresponding to the 
various modes involved will all have waves propagating without attenua- 
tion: the introduction of lines corresponding to evanescent modes allows 
another class of problem to be solved. The propagation characteristics of 
artificial dielectrics (3) or corrugated waveguides are often calculated by 
means of an equivalent circuit: this method is usually limited to the 
situation when the reactive fields associated with the elements of the 
medium do not interact. The introduction of transmission lines corre- 
sponding to the most important of the evanescent modes enables this 
limitation to be removed. 

In this paper some general results applicable to systems of coupled 

(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 2, 1958] 
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transmission lines are presented. It is shown that when the lines are 
coupled at regular intervals by identical coupling networks a set of modes, 
each capable of a separate existence, can be found. The properties of these 
modes show certain similarities to those arising in conventional waveguides. 

We consider a system consisting of N transmission lines 1, 2,..., N, 
connected together by networks X spaced at equal intervals a along the 
lines as in Fig. 1. The propagation constants of the lines are y,, y9,.--. ¥x 
and their characteristic impedances are Z,, Zg,..., Z,, respectively. We 





< SS ee 


Fia. 1 


restrict the analysis to loss-free lines so that y is either purely imaginary 
for a line on which propagation can occur or purely real for evanescent 
lines. Similarly, Z is positive real when propagation can occur and is 
imaginary when the line isevanescent. The sign of the imaginary impedance 
for an evanescent line is positive if the infinite line stores magnetic energy 
and is negative if it stores electric energy. The networks X are assumed 
to be purely reactive. This restriction to purely loss-free systems does not 
exclude any of the problems mentioned in the introductory paragraph: 
the extension of the analysis given here to the general case when loss is 
present in the lines, the coupling networks, or both, complicates the algebraic 
treatment but does not require any additional principles. 

The network of Fig. 1 can be split into a number of identical sections such 
as that shown in Fig. 2. The method of analysis is formally the same as 
that used in filter theory, namely, that we seek solutions which have the 
property that the voltages and currents at the output terminals B differ 
from those at the input terminals A only by a common factor which has 
the significance of a propagation constant per section. A preliminary step 
is to obtain the general relation which must be satisfied between the output 
and input voltages and currents. This can be carried out by using transfer 
matrices. 
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2. The transfer matrix for the basic section 

The voltage at any point in Fig. 2 is denoted by V, p, where n is the 
number of the line across which the voltage is developed and P is the 
position at which it is observed. Similarly, J, » is the current at position P 
on line n, the direction of positive current flow being always from left to 
right. 





Fia. 2 


Between A and C propagation on any line occurs independently of all 
the others and from transmission line theory we have that 
V4 = Vi, cosh by, 4+ -Z, Ic sinh hy, a (1) 
TA —_ (V,, «sinh by, @)/Z,,+-1,,¢ cosh ty, @ ; 
Equations similar to (1) apply for each of the N lines, and the results may 
be collected in matrix form to give 
(Vs) _ C ZS]{V\ 
\1,J ZS C |\I,J’ 
where V, is the column matrix {V,4,V24,...,Vy4} with L,, Vo, and I, 
similarly defined, and C, Z, S are diagonal matrices in which the diagonal 


(2) 


elements are 

cosh }y,a, cosh}y,a, ..., cosh}yya 

By Be. i Be 

and sinh }y,a, sinh}y,a, ..., sinh}yya 
respectively. We may note that since diagonal matrices commute, the 
order of Z and S in the product is immaterial, and also that Z~' will always 
exist. Further, all the elements of C are real since each y is either purely 
real or purely imaginary and all elements of ZS and @-'S are imaginary 
because of the restrictions imposed on y, and Z,, by virtue of the loss-free 
nature of each line. 
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The connexion between D and B is identical to that between A and C 


so that r 
(Vp\ _ C ns bag (3) 
\I,) Z‘'S C|\I;,!) 
where the column matrices V,,, etc., are defined in a similar way to V,. 
Terminal positions C and D are connected by the reactive network X, 
so that ' 
(Vol — jx{ Fe | (4) 
\V,J | —I,) 
where X is a real symmetric matrix with 2N rows and columns. The minus 
sign associated with I, appears since the form of (4) requires that the 


positive direction of current flow be towards the network. Equation (4) can 
be written in partitioned form as 


Ne) fA BY Be) 6 
\Vp) B D \—1!’ 
A, B, D being real matrices with V rows and columns. A and D are each 
symmetrical. 
A further restriction will be imposed in that X is assumed to be sym- 
metrical with respect to terminals C and D, which requires that 


A=D; B = B. (6) 
Equation (5) is now rewritten to relate V, and I, to V, and I,: 
V. = jAI,—jBI),, (7) 
V, = jBIp—jAlp. (8) 
Manipulation of these equations by matrix algebra gives 
V, = AB-"V,+j[AB“A—BII,, (9) 
I, = —jB-"V,+B-Al),. (10) 


Since B is an impedance matrix, the reciprocal B-' exists. The combina- 
tion of (9) and (10) gives 


(Vo\_ [AB j(AB-'A—B)]/V,\ (11) 
Io} iB BA it }’ 


and from (2), (3), and (11) we have finally 


= (55 Se LS Ae). an 
‘in lss cla es lize char ® 


which relate the voltages and currents at the input and output terminals, 
This can be written as 


D 


as = 
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the elements of the transfer matrix being found by expanding the matrix 
product : 

T = CAB-'C+SZAB-'SZ-!—jSZB-'C+jC(AB-'A—B)SZ-", 

U = CAB-'SZ+SZBAC—jSZB-'SZ+jC(AB-'A—B)C, (14) 

V = SZ"*AB-'C+ CB-'ASZ—jCB“'C+jSZ-(AB-'A— B)SZ-". 
Since A, B, C, S, and Z are all symmetrical, it is easily verified that U 
and V are symmetrical. Further, from known properties of the matrices 
on the right of the equations, it may be shown that T is real while U and 
V are purely imaginary. The following relations between T, U, and V can 
also be verified by straightforward matrix multiplication 


TU = UT, (15) 
VT = TV, (16) 
UV = T?—I. (17) 


3. Modes of propagation 


A single mode of propagation can be defined as occurring when the 
voltages and currents at B (Fig. 2) differ from those at A by some constant 
factor which can be written exp(—Ta). Let a, B be the corresponding 
voltage and current matrices at A, so that 


V, = exp(—Ta)a; I, = exp(—Ta)B. (18) 
Substitution in (13) gives 
s}—ew-rof gh 
which becomes on expansion 
a = exp(—TIa)(Ta+ UB), (20) 
® = exp(—a)(Va+ Tp). (21) 
From (20) and (21) we get 
UB = (exp(l'a)I—T)a, (22) 
(I—exp(—Ta)T)p = exp(—Ta)Va. (23) 
Pre-multiplying (23) by U and using (15) and (17) we get 
(I—exp(—Ta)T)UB = exp(—Ta)(T?—I)a. (24) 


Substituting from (22) for UB and simplifying we find 
Ta = cosh(Ta)a (25) 
from which we see that non-trivial solutions arise only if cosh([a) is a 


latent root of T. There are N such roots giving N possible modes of 
propagation for the system of N coupled transmission lines. 
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It will now be shown that all the latent roots of T are real; this is clearly 
the case using the standard result that the latent roots of a real symmetric 
matrix are real. We select a real matrix M such that 


MM = jv, (26) 

which is always possible since jV is a real symmetric matrix. 
Using (16) we have MMT = TMM, (27) 
whence MTM-! = MT, (28) 


so that MTM~' is a real symmetric matrix. Now let 
a = M-ly (29) 
and substitute in (25) giving 
TM~'y = cosh(Ta)M-'y, 
i.e. MTM-'y = cosh(Ta)y. (30) 


This shows that cosh(Ta) is a latent root of a real symmetric matrix and 
is therefore always real. The effective propagation constant [’ must there- 
fore be purely real or purely imaginary, showing that the mode either 
propagates in a similar way to a normal waveguide mode or is evanescent. 

Waveguide modes are typified by the distribution of electric and mag- 
netic fields in a plane transverse to the direction of propagation. Corre- 
spondingly, the modes here will depend on the distribution of voltages on 
the transmission lines at planes such as A or B (Fig. 2), i.e. on the column 
vector a which is the latent vector of T corresponding to the latent root 
cosh(T'a). We may number the modes, a convenient method being to order 
them in terms of increasing values of the latent roots. Then, a, is the 
latent vector corresponding to the latent root cosh(I,a), i running from 
1 to N. The elements of a; can be determined, apart from an arbitrary 
constant, by solving the linear homogeneous set of equations given by (25). 
Since T is a real matrix and cosh(Ija) is real, the constant can always 
be selected to make all the elements of a; real. This implies that when 
a single mode is propagating along the system of transmission lines, the 
voltages at any plane such as A or B are all in phase. 

We have tacitly assumed above that the latent roots are all different, 
but in degenerate cases two or more roots may be identical. The latent 
vectors are then no longer uniquely defined. To preserve a unified treatment 
we suppose that the vectors are selected so that the complete set of latent 
vectors remains an orthogonal one: this can always be done. The rest of 
the analysis is unaltered. Degenerate modes in waveguides are handled in 
essentially the same way in that the selection of mode patterns corresponding 
to the degenerate modes is made to maintain orthogonality. 

We now require the corresponding current vector which is given by 
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(22) or (23). The most convenient form for B; is obtained by repeating the 
analysis of (20) to (25) to eliminate a; when it is found that 


| TB, = cosh(I,a)B,. (31) 
If this is substituted in (23) and the equation simplified, we find 
8B, = cosech(T,a)Va; (32) 


so that 8, is completely determined. For a propagating mode, [T, and 
cosech(I,a) are imaginary, so that all the elements of 8, are real, since V 
is purely imaginary and a; is real. Thus, for such a mode the currents and 
voltages on all the lines are in phase at planes such as A or B. When the 
mode is evanescent, a similar argument shows that all the elements of B; 
are imaginary, the currents then being in phase quadrature with the 
corresponding voltages. It may be observed that similar results apply to 
the transverse distributions of electric and magnetic field strengths in 
waveguide modes. 


4. Orthogonality properties of the modes 
The essential orthogonality property satisfied by waveguide modes is 
that when several modes are propagating simultaneously the total power 
flow equals the sum of the power flows for the modes taken one at a time. 
A similar result holds in the present case. We start by obtaining an 
expression for the power flow in the most general case when all the modes 
are excited. We note first that the modes can travel from right to left as 
well as from left to right, the only change needed being in the sign of Ij. 
Suppose then that the amplitudes and phases of the components of mode i, 
as measured at plane A, are given by the complex numbers P; and Q;,, 
respectively. This implies that: 
V, = (P.+- Qa, (33) 
1, = (P—Q)Be (34) 
the reversal of sign in the second term of (34) arising from the occurrence 
of cosech(I,a) in the expression for B; (equation (32)). The power flow 
towards the right at plane A is the sum of the power flows into each of 
the individual lines, i.e. 
power flow = bre Via It, = bre WT, (35) 
n=1 


~ 


V, being the row vector obtained by transposing the column vector V,, 
and the asterisk denoting the complex conjugate. 
From the expressions in (33) and (34) coupled with (32), we have, 
power flow = }re(P;+ Q,)( P¥ — Qf)a; V*a, cosech(I'¥ a). (36) 


The corresponding expression for the power flow on a transmission line 
5092 .42 F R 
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of characteristic admittance »;, when there are forward and backward 
waves given by the complex numbers P, and Q,, respectively, is 
power flow = }re(P;+ Q,)( P?— Q?)17. (37) 


Comparison of equations (36) and (37) shows that the effective admittance 


for mode i is 1 = &, Va, cosech(T;a). (38) 


From the real and imaginary properties of the various quantities on the 
right, it is found that 7, is real or imaginary according as the mode is 
propagating or evanescent respectively. Further, since the elements of a, 
involve a common arbitrary constant, the values of 7; are arbitrary apart 
from their real or imaginary character. Again, we may note a similarity 
to waveguide modes for which there is no unique definition of impedance. 
When all modes exist simultaneously (33) and (34) are generalized in 
that the right-hand sides now involve summations with i ranging from 1 
to N. These equations may be written most concisely by defining square 
matrices R and S in which the ith columns are a, and B,, respectively. 
The expressions for the voltages and currents at A are then 
V, = RP+RQ, (39) 
I, = SP—SQ, (40) 
P and Q being column matrices formed from the elements P,; and Q; 
respectively. With the above definitions for R and S, (32) can be extended 
to give S = VRA. (41) 
A being a diagonal matrix, in which the diagonal elements are cosech(T, a) 
(¢ <= ], 2...., N). 
The expression for the power flow when all N modes exist is the generali- 
zation of (36) 
power flow = }re[P +Q|RV*RA*(P* —Q*}. (42) 
The orthogonal property is demonstrated by showing that RV*R is a 
diagonal matrix: to do this we use the expression for «; in terms of y; given 
in (29). Again, this can be extended to give 
R = M~y. (43) 
The columns of y are the latent vectors of the real symmetric matrix 
MTM-~' so that by a standard result 
¥y = a diagonal matrix, (44) 
giving RMMR = jRVR = a diagonal matrix. 
The values of the diagonal elements follow from (38), whence 
RVRA = H, 


H being a diagonal matrix with diagonal elerents », (¢ = 1, 2...., 
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Substitution in (42) gives at once 


N 


power flow = }re ¥ (P+ Q,)nf(P7—Q?), (47) 


=1 
showing that the total power flow is obtained by adding the contributions 
of the modes, considered one at a time. 
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SUMMARY 


The propagation of electromagnetic waves in an inhomogeneously filled wave- 
guide is discussed and it is shown that, when the phase velocity of the wave down 
the guide is equal to the velocity of propagation in one of the dielectric media, then 
the longitudinal component of magnetic field is constant across that dielectric. 
Limits are also given for the propagation constants of TE and TM waves. 


1. Introduction 


Ir has been stated (1) that it is possible under certain conditions to propa- 


gate TEM waves down waveguides which are filled with more than one 
dielectric, the phase velocity being the velocity of light in one of the com- 
ponent dielectric portions of the guide. This statement has been quoted 
in a summary of the work which had then been done on inhomogeneously 
filled guides (2). 

It has been pointed out (3) that the existence of a wave having the 
velocity of light in the medium does not of necessity imply vanishing of 
the longitudinal components of field strength, but only of its transverse 
derivatives. That is to say, propagation with the velocity of light in the 
medium does not always involve a TEM wave. In general, propagation 
of this character will involve finite Z, and H, (z being the longitudinal 
direction), but H, and E, will be constant across that portion of the guide 
such that the velocity of light is the velocity of propagation. Under 
certain conditions, however, the quantities Z, and H, will vanish, and 
TEM waves will occur. When a field is such that H, does not have trans- 
verse variation, within a region it will be said to be constant longitudinal 
magnetic (CLM). Similarly, if 2, does not vary, the field is said to be 
constant longitudinal electric (CLE). (It will be noted that CLE fields 
reduce to TE fields when the region is contiguous to a metal wall (Z, = 0).) 
Conditions for the existence of CLM, CLE, and TEM waves in a certain 
rectangular guide will now be discussed. 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 2, 1958) 
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2. Bounds on propagation constant in an inhomogeneously filled 
waveguide 
It may be shown that if a TE or TM wave is propagating down a wave- 
guide, the electromagnetic constants of which vary over its cross-section, 
the value of 8”, the square of the propagation constant, must be less than 
the maximum value of k? (= we) across the cross-section. The proof runs 


as follows. For a waveguide whose dielectric properties vary across its 
cross-section (4) 


{ {uJH|*—<|E}dW = 0, (2.1) 


Ww 


where W is the waveguide cross-section, » is the magnetic permeability, 
and «¢ the dielectric constant. 
An arbitrary TE field can be expressed in the following form : 


E = —wwA, (2.2) 
H = |¢y,+ipk}, A, (2.3) 
ub 
where A,= 0. Y, is the transverse gradient operator and a factor 


exp{i(wt+-8z)} is assumed. k is unit vector in the z-direction. It follows 
from equations (2.2) and (2.3) that 


E — iwA, (2.4) 
fi — ; (Vv, —iBk) A, (2.5) 


where the tilde denotes the complex conjugate value. Substituting into 
equation (2.1) it follows that 


[| 00. + 48k)» A}.{00, —iBk) A B}—eut.A dW=0. (2.6) 


Ww 


This simplifies to 


[ [Ei AAP+BIAI}—eutiAl| aw = 0, (2.7) 
W 
[ “tw. AA|?+(f2—k)|A)?] dW = 0. 

pe 


. 


Ww 


Now if 8? > k® always, the integrand is positive definite, and equation 
(2.7) can only be satisfied by the trivial solution of zero A. Thus it follows 
that 6? must be less than the maximum value of k?. The proof for a TM 
wave follows with very slight alterations. 
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3. Rectangular guides 

Consider the waveguides whose physical characteristics are of the form 
indicated in Fig. 1. For simplicity it will be assumed that the system is 
symmetric about z = 0. Clearly the waves such that E, is antisymmetric 
about x = 0 will be the same as those in a waveguide with a metallic 
wall at x = 0. 

The fields which exist are assumed to be of the form F(x)exp{i(wt+8z)}. 
It is the object of this investigation to find the condition for fields such 


w—a-l —>e— 21 —e—a-l > 
Eis fy €2> #2 ee a 




















Fie. l 


that H, is constant within one of the regions 0 < |x| <l orl < |x| <a. 
Obviously a wave with EZ, constant would not be possible because E£, 


<= 


vanishes on the walls. A field which is independent of y is given by 
@E, 


Gat t+ (uew*—B2)E, = 0, (3.1) 


uv 


i ee (3.2) 


z 


cope 

eee. (3.3) 
iw dx 

The boundary conditions which must be obeyed are 


E, H,=9 at |z| =a, (3.4) 


| 
and E,, wH,, H, are continuous at |x| = 1. It is only necessary to consider 
the fields in 0 < 2 < a, as the fields in —a < x < 0 can be dealt with 
by using symmetry or antisymmetry. It will be noticed that this wave 
is of TE character. 


Case (a). B® = wp, €, 
By section 2, it follows that a wave is only obtained when 


He € > fy &- 
Ini <2 <a, 
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Thus Z, is linear, and in order that (3.4) may be satisfied 
E,=E,(l—z/a) in l<2z<a, 
Ey 


whence H =- : 
twp, a 


IndO<2x<l, OO atna~mioe = 0, 
dx? . . 
the solution of which is 
E, = Asin px+ Becos pz, 
—twp, H, = Apcos px— Bpsin pz, 
where p? = w*(py €g—py €&)- 


If H, is antisymmetric, the matching conditions at x = I give 
B{1—*) = Beos pl, (3.13) 


Me, PB ings (3.14) 
Wwp,4 Wh» 


and, eliminating EZ, and B, it follows that 


pltan pl = 4 : ) (3.15) 


py \a—l 
and this gives the frequencies at which waves of the character considered 
are possible. 


If on the other hand H, is symmetric, the matching conditions are 


B{1—) = Asin pl, (3.16) 


8. = ~ cos pl. (3.17) 
twp, a Who 
Eliminating 2, and A, it follows that the equation determining the fre- 
quency is given by I 
—plcot pl = Ha } (3.18) 
My 


a—l 
Case (b). B® = wp, €, 

By section 2, it follows that a wave is only obtained when p,«€, < p, €. 
In 0 < x < 1, the fields are given by 
(3.19) 


(3.20) 


whence (3.21) 
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d?E 
Vi .2 avs j 
az? + w*(py €; Hy €)E, = 0, 


The solution of (3.23) is given by 

E, = Acosqx+ Bsinqz, 
where q? = wy €;— pe €). 3.26) 
To satisfy the boundary condition at x = a, we must have 


E, = Csing(a—2), (3.27) 


H, = - Cq cos q(a—2). (3.28) 
twp 
Now if H, is antisymmetric, A vanishes, and it follows that H, vanishes; 
and that in this particular case the wave is TEM in 0 < x < l. 
The frequency equation is obtained by the application of the matching 
condition at z = / for H,, and the frequencies are given by 


0 = cosq(a—l), (3.29) 

whence q(a—l) = (2n—1)}x,_ n integral. (3.30) 
Consider now the case of symmetry in H, about x = 0. When this 
happens Z,, vanishes at x = 0, and the system is the same as if a metallic 
wall were inserted at x = 0. It follows that the problem is the same as 
that discussed in equations (3.16) to (3.18) with the media 1 and 2 inter- 


changed and the lengths /, a—/ also interchanged. Thus it is sufficient to 


consider only the case baile > Pits (3.31) 


4. Numerical results 
In figs. (8. 1-2) and (8. 1-3) of (5) a number of graphs have been drawn 


giving the propagation wavelengths for inhomogeneously filled guides of 
the character indicated for the case 


Be = Hy = Ho» (4.1) 
€, = 2:45€, = 2-45¢,. (4.2) 


In these diagrams, the fields which are of interest to this paper are those 
for which d,/A, = 1. (4.3) 
Fig. (8. 1-2) gives the case of antisymmetric H, and in this figure a is the 
complete guide width. Fig. (8. 1-3) gives the case of symmetric H, and is 
for a guide of width a with a dielectric slab at one side of it. From 
section 3, this also gives the symmetric H, modes in the guide of width 2a. 
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In order to get a general idea of the manner in which the fields behave, 
the case a = 2/1 will now be considered in detail. The electromagnetic 
constants will be given by (4.2) and thus 


p = 1-2047k, where k2 = wy €. (4.4) 
The frequency equations reduce to 


xtany = 1 (4.5) 


for antisymmetric H,, and ycoty = —1 (4.6) 


for symmetric H,, where pl = x. (4.7) 


The solutions of equation (4.5) and equation (4.6) interlace one another 
as might be expected. 
The solutions of equation (4.5) are 
x = 08604, 3-4256, 6-4373.... (4.8) 
and of equation (4.6) 
x = 2-0288, 4-9131, 7-9737,.... (4.9) 
It follows from equation (4.7) that the free space wavelength A, corre- 
sponding to the fields is given by 
Ay = (U/x)1+2047. 27 
= 7-570(I/,). (4.10) 
The first antisymmetric mode is therefore given by A,/l = 8-798 and 
the first symmetric mode by A,/l = 3-731. 
The behaviour of H, may be obtained from equations (3.9) and (3.12). 
For antisymmetric H,, 
H,=sinpxr (0 < |z| <1) 
sinpl (l< |x| <a). (4.11) 
For symmetric H,, 
H,=cospx (0 < |x| <1) 


cospl (1 < |x| <a). (4.12) 


Tables for the behaviour of H, in the first two modes follow for the case 
@ = 2. 





x/l ° orl o°2 o-3 o-4 o'5 06 o'7 o'8 o'9 10 





ist mode : ° 0086 | O°171 | O'255 | 0°337 | 0°417 | 0°494 | 0566 0°635 0699 0758 
2nd mode. I 0980 | 0-919 | 0°820 | 0°688 | 0-528 | 0-346 | o-150 | —o°052 | —0°252 | —0-442 






































A graph of the behaviour of H, across the guide for the first two CLM type 
modes is given in Fig. 2. 
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5. Discussion 

It is of interest to discuss how the CLM waves are related with other 
types of waves which may occur in a waveguide. CLM waves are in fact 
particular cases of Longitudinal section electric (LSE) and Longitudinal 
section magnetic (LSM) (6) waves in which there is no variation in the 
y-direction, and the CLM wave is in fact a wave which is LSE, LSM, and 
TE simultaneously, this being a condition which can only be satisfied 
when there is no variation in the y-direction, and when the velocity of 
propagation of the wave along the guide is the same as the phase velocity 
of light in one of the component dielectric portions. 

Consider now a guide 0 <4 <a, 0<y <b, such that the electro- 
magnetic constants in region p (x,_, << 4% << %);% = 0, x, = a) aree,, py. 

A field in region p which is independent of y is given by 


a*E. ee 
att ine o® BE, = 0, 
ee Sh. 
iw, dx 
Equation (5.1) may be rewritten 
dE, 
dx? 


— BL E, = 9, 


a: ae ies 
where f% is not necessarily positive. 
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The solution of equation (5.3) is given by 
E, = A, cosh B,, x+ B,(sinh 8, x)/B,,, (5.4) 
where (sinhf,,x)/B, =x if B, = 9, (5.5) 
—iwH, = pn, {A, , sinh 8, x+ B, cosh B, 2}. (5.6) 
The boundary conditions are that Z, and H, are continuous across z,, 
and that Z, = 0 onz = Oand z= a. 
Hence A, = 0, from the boundary condition at x = 0, and 
A,, cosh 8, a+ B, (sinh 8,,a)/B, = 0 (5.7) 


from the boundary condition at x = a. 
The boundary conditions at x, (1 < p < n—1) give 


A, cosh B,, x, + B,(sinh B, x,)/B, 
= A,,, coshB,.,2,+B,.,(sinhB,.,%,)/Bpi1, (5.8) 


up (A,B, sinh B, x,+B, cosh B, z,) 

= ppi(A,B8,.,8inhB,,,27,+B,,,coshB,,,2,). (5.9) 
Now if the value of 8 is such that 8, vanishes, then, provided that the 

conditions implied by section 2 apply, a wave exists such that 
E, = A,+ B,2, (5.10) 
—iwpu, H, = B,, (5.11) 
and in this case the CLM wave exists. In particular, if B,, vanishes the wave 
will be TM, and since it is already TE, is therefore a TEM wave. It will 
be seen that the equations (5.7) and (5.8) between them determine the ratios 
of (2n—1) constants. It follows that when f is such that 8, vanishes, 
there is always a CLM wave, provided that the conditions of section 2 
apply, but that the existence of a TEM wave depends on B,, vanishing, 
a contingency which clearly happens only under certain critical relations 
between the x,, ,, €,. An example of this has been discussed in section 3. 
It will be seen that similar phenomena occur in circular waveguides 
which either have a central core of dielectric, or are lined with dielectric 
on the inside of the bounding metal walls. In the first case, it will be 
possible to propagate a TE wave which is CLM in the region between 
core and guide walls, and in the second case a TE wave which is CLM in 
the central portion of the guide. In addition, in this latter case there is 
no metal in the central portion of the guide and so it is not necessary for 
E, to vanish therein and a TM mode which is CLE propagates at the 
appropriate frequency. Thus if a field is required such that a component 
of magnetic field does not vary in the plane transverse to it, but varies 
with a particular frequency, it is possible to modify suitably either a 
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circular or rectangular guide with dielectrics. If on the other hand the 
same requirement be imposed on an electric field component, it is neces- 
sary to use a circular waveguide. 
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PROPAGATION IN A GYRATIONAL MEDIUM 
(ADDENDUM) 


By LL. G. CHAMBERS (University College of N. Wales, Bangor) 


[Received 20 May 1957] 


SUMMARY 


A previous paper (1) by the author discussed the electromagnetic properties of 
gyrational media, but reserved the question of ideally gyrational media for further 
discussion. The present paper provides the further discussion required, and the 


enumeration of sections, equations, and reference follows on consecutively that of 
the previous paper. 


9. Ideally and perfectly gyrational media 
WHEN pe—EF = O (9.1) 
the medium is ideally gyrational, and the constitutive equations are of the 


form D = cE+a,(pe)H, (9.2a) 
B = a 1\/(we)E+puH. (9.2 b) 
Clearly for such media, k? (1, p. 361) vanishes anc 
D = ay(e/u)B. (9.3) 
It also follows that it is impossible to solve the constitutive equations 
(9.2) for E and H. (It will be noted that if «? is unity, the quantity 
7 = tw(C—E€) which was introduced in section 2, and which gives rise to 
the doubly refracting property in the non ideally gyrational case vanishes, 
and the medium may be said to be perfectly gyrational.) Utilizing equation 
(9.3), Maxwell’s equations become 
VAH—iwa/(e/u)B = 0, (9.4a) 
VAE+ioB = 0. (9.4b) 
From these equations it follows that 
V/A {pH +071 (ue) E}+-iwfat, (we) —ay/(we)}B = 0 
and, using the constitutive equations, this may be rewritten as 
VA B+ 7B = 0. 
Also, from equation (9.4b), V.B=0. 
Taking the curl of equation (9.5) we have 
VA(VAB)+7V AB = 0, 
whence V(V.B)—V*B—7*B = 0, 
and using equation (9.6), 
V°B+-7°B = 0. 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 2, 1958) 
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It follows, therefore, that B propagates with speed (w/) and there is, in 
this case no double refraction. Furthermore it follows from equation 
(9.5) that the waves are circularly polarized in a definite sense. If» vanishes, 
the distribution of B and D is identical with a static distribution. 

It is convenient to attempt to express the field components in terms of 
scalar and vector potentials. Firstly, it is clear from equation (9.5) that if 
B = VAA, then the vector potential is given by 


A —7 1B (9.8) 
and A satisfies the same differential equation as B. Also, from equations 
(9.4) and (9.5), VA {nH +iway(e/p)B} = 0, 

VA (nE—iwB) = 0. 


It follows that H = —iwan-'(e/n)B—VU, (9.9 a) 


E = iwn'B—VV, (9.9b) 


where U and V are arbitrary point functions, though related by equation 
(9.2 b). It follows that 
U —aW(e/p)V. (9.10) 


The equation (9.9 b) is in fact, using equation (9.8), the well-known relation 
E = —iwA—VV. (9.1la) 
The remaining field components are 


H = iway(e/u)A+aV(e/n)VV, (9.11 b) 


D = —ay(e/u)nA = av(e/p)WAA, (9.11 ¢) 
B=—,nA=VAA. | (9.11 d) 


A satisfies the differential equation (9.11d) and V is arbitrary. The 
arbitrariness of V corresponds to the impossibility of solving the consti- 
tutive equations for E and H. 

The expression of field components in terms of potentials indicated above 
depends on the non-vanishing of 7, and hence this case (which corresponds 
to a? unity) still remains to be treated. Clearly from equation (9.4 b) 


B= VAA, (9.12 a) 
E = —iwA—VV. (9.12 b) 
From equation (9.3) D = av(e/u)WAA, (9.12 ¢) 
and from equation (9.2 b) 
H = p B—aY(e/p)E 
BOW A A+iwa)(e/p)A+a/(e/n)0V (9.12d) 
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using the fact that a* is unity. Substituting equation (9.12 a) in equation 
(9.5), it follows that A satisfies the differential equation 


VA(VAA) = 0; (9.13) 
V is again arbitrary. 

The solution of equation (9.7) presents no difficulty whether » vanishes 
or not, being in fact a standard equation of electromagnetism. Accordingly, 
it is not deemed necessary to indicate the different types of field which can 
arise. 
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CORRIGENDA 


J. S. LOWNDES: ‘A TRANSIENT MAGNETIC DIPOLE SOURCE 
ABOVE A TWO-LAYER EARTH’ 


(Vol. X, Pt. 1, 1957) 


. 80, the second of equations (2): ‘H;,’ should be ‘H;,’ 
equation (23): ‘M’ should be ‘2M’ 


) should be ‘2(2) 


Fig. 1 now gives the values of A(t)+ (S)ere’(S,) 


5, the expression for A(t): (=, 


‘ \3? ‘ 3° 
}, the expression for B(t): s(25) should be (=) 
Fig. 2 now gives the values of B()—2() ert’(=) 
ites a 2¢4/ 2¢4 
. 88, the second of equations (43): ‘ F(t)’ should be ‘—¢t~! F(t)’ 
. 89, the expressions for D(t): ‘(2az*)-?’ should be ‘(2z2z*)!” 
‘ } ,’ ‘ , 
. 89, equations (44) to (47): —= should be +3 
. 89, equation (47): ‘MM’ should be *‘ Mit-!” 


The author regrets that the above corrections were not made before the publica- 
tion of the paper. 
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